Changeset 4071 for trunk/LMDZ.COMMON/libf/evolution/pem.F90

Timestamp:

Feb 16, 2026, 10:28:56 AM (12 days ago)

Author:

jbclement

Message:

PEM:

• Making the computation of ice tendencies more reliable by doing it after 'read_startpem' to know exactly where perennial ice is (no matter if there is a "startpem.nc" or not). Moreover, when there is no "startpem.nc", location of perennial ice depends now on 'watercaptag' and on huge amount of frost. This prevents negative ice tendency while there is no ice which can happen with weird PCM inputs (i.e. 'watercaptag' = F & 'watercap' /= 0 & no "stratpem.nc").
• Few small safeguards throughout the code.

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/pem.F90 (modified) (7 diffs)

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

@@ -48 +48 @@
 use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs, build4PCM_perice
 use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice, build4PCM_tsurf
-use tendencies,         only: evolve_tend_co2, evolve_tend_h2o
+use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_tend_h2oice
 use tracers,            only: adapt_tracers2pressure, build4PCM_tracers, nq
 use utility,            only: real2str
@@ -77 +77 @@
 real(dp)                              :: preff4PCM       ! Reference pressure [Pa]
 ! Ice-related:
-real(dp),    dimension(:,:), allocatable :: h2o_ice                    ! H2O ice [kg.m-2]
-real(dp),    dimension(:,:), allocatable :: co2_ice                    ! CO2 ice [kg.m-2]
-real(dp)                                 :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2]
-real(dp)                                 :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2]
-logical(k4), dimension(:,:), allocatable :: is_h2oice_ini              ! Initial location of H2O ice
-logical(k4), dimension(:,:), allocatable :: is_co2ice_ini              ! Initial location of CO2 ice
-logical(k4), dimension(:,:), allocatable :: is_co2ice_flow             ! Flag for location of CO2 glacier flow
-logical(k4), dimension(:,:), allocatable :: is_h2oice_flow             ! Flag for location of H2O glacier flow
-real(dp),    dimension(:,:), allocatable :: h2o_ice4PCM                ! H2O ice reconstruction to feed back into PCM [kg.m-2]
-real(dp),    dimension(:,:), allocatable :: co2_ice4PCM                ! CO2 ice reconstruction to feed back into PCM [kg.m-2]
-logical(k4), dimension(:),   allocatable :: is_h2o_perice              ! Location of H2O infinite reservoir, called 'watercaptag' in PCM
-logical(k4), dimension(:,:), allocatable :: is_co2ice_disappeared      ! Flag to check if CO2 ice disappeared at the previous timestep
+real(dp),    dimension(:,:),   allocatable :: h2o_ice                    ! H2O ice [kg.m-2]
+real(dp),    dimension(:,:),   allocatable :: co2_ice                    ! CO2 ice [kg.m-2]
+real(dp)                                   :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2]
+real(dp)                                   :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2]
+logical(k4), dimension(:,:),   allocatable :: is_h2oice_ini              ! Initial location of H2O ice
+logical(k4), dimension(:,:),   allocatable :: is_co2ice_ini              ! Initial location of CO2 ice
+logical(k4), dimension(:,:),   allocatable :: is_co2ice_flow             ! Flag for location of CO2 glacier flow
+logical(k4), dimension(:,:),   allocatable :: is_h2oice_flow             ! Flag for location of H2O glacier flow
+real(dp),    dimension(:,:,:), allocatable :: minPCM_h2operice              ! Minimum of H2O perennial ice over the last PCM year [kg.m-2]
+real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice              ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2]
+real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost               ! Minimum of H2O frost over the last PCM year [kg.m-2]
+real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost               ! Minimum of CO2 frost over the last PCM year [kg.m-2]
+real(dp),    dimension(:,:),   allocatable :: h2o_ice4PCM                ! H2O ice reconstruction to feed back into PCM [kg.m-2]
+real(dp),    dimension(:,:),   allocatable :: co2_ice4PCM                ! CO2 ice reconstruction to feed back into PCM [kg.m-2]
+logical(k4), dimension(:),     allocatable :: is_h2o_perice              ! Location of H2O infinite reservoir, called 'watercaptag' in PCM
+logical(k4), dimension(:,:),   allocatable :: is_co2ice_disappeared      ! Flag to check if CO2 ice disappeared at the previous timestep
 ! Surface-related:
 real(dp), dimension(:,:), allocatable :: tsurf_avg           ! Average surface temperature [K]
@@ -179 +183 @@
 allocate(tsoil_avg(ngrid,nsoil,nslope),tsoil_ts(ngrid,nsoil,nslope,nday),h2o_soildensity_avg(ngrid,nsoil,nslope))
 allocate(q_h2o_ts(ngrid,nday),q_co2_ts(ngrid,nday))
-allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope))
+allocate(minPCM_h2operice(ngrid,nslope,2),minPCM_co2perice(ngrid,nslope,2),minPCM_h2ofrost(ngrid,nslope,2),minPCM_co2frost(ngrid,nslope,2))
 
 call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, &
-                    q_h2o_ts,q_co2_ts,d_h2oice,d_co2ice)
+                    q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
 
 ! Initiate soil settings and TI
@@ -207 +211 @@
                    h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
 deallocate(tsurf_avg_yr1)
+
+! Compute ice tendencies from yearly minima
+allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope))
+call print_msg('> Computing surface ice tendencies')
+call compute_tendice(minPCM_h2operice + minPCM_h2ofrost,h2o_ice(:,:) > 0._dp,d_h2oice)
+call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)))
+call compute_tendice(minPCM_co2perice + minPCM_co2frost,co2_ice(:,:) > 0._dp,d_co2ice)
+call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)))
+deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
 
 ! Save initial set-up useful for the next computations
@@ -290 +303 @@
         ! Conversion to surface ice
         call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
+        ! Balance H2O ice reservoirs
+        allocate(d_h2oice_new(ngrid,nslope))
+        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)
+        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)
+        deallocate(d_h2oice_new)
     else
         zlag(:,:) = 0._dp
         call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
         call evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
-    end if
-
-    if (do_layering) then
-        allocate(d_h2oice_new(ngrid,nslope))
-        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)
-        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)
-        deallocate(d_h2oice_new)
     end if
 
@@ -406 +417 @@
 
     ! Evolve the tendencies
-    call evolve_tend_co2(d_co2ice_ini,co2_ice,emissivity_PCM,q_co2_ts_ini,q_co2_ts,ps_ts,ps_avg_glob_ini,ps_avg_glob,d_co2ice)
-!~     call print_msg("> Updating the H2O sub-surface ice tendency due to lag layer")
+    call evolve_tend_co2ice(d_co2ice_ini,co2_ice,emissivity_PCM,q_co2_ts_ini,q_co2_ts,ps_ts,ps_avg_glob_ini,ps_avg_glob,d_co2ice)
 !~     if (do_layering) then
 !~         do i = 1,ngrid
 !~             do islope = 1,nslope
-!~                 if (is_h2oice_sublim_ini(i,islope) .and. h2oice_depth(i,islope) > 0._dp) call evolve_tend_h2o(h2oice_depth_old(i,islope),h2oice_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
+!~                 if (is_h2oice_sublim_ini(i,islope) .and. h2oice_depth(i,islope) > 0._dp) call evolve_tend_h2oice(h2oice_depth_old(i,islope),h2oice_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
 !~             end do
 !~         end do
@@ -417 +427 @@
 !        do i = 1,ngrid
 !            do islope = 1,nslope
-!                call evolve_tend_h2o(icetable_depth_old(i,islope),icetable_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
+!                call evolve_tend_h2oice(icetable_depth_old(i,islope),icetable_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope))
 !            end do
 !        end do