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

Timestamp:

Apr 3, 2026, 4:34:51 PM (8 days ago)

Author:

jbclement

Message:

PEM:

• Deletion of 'flux_ssice' ('zqdsdif_tot') from the "startfi.nc" as it is not needed.
• Using the yearly average flux for the sublimating subsurface ice instead of the value at last timestep.
• Keeping a clear separation between the subsurface ice flux and the surface ice tendency.
• Making sure that subsurface ice depth is well given to the PCM at the end of the PEM (ice table VS layering).

JBC

File:

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

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

@@ -52 +52 @@
 use surf_ice,           only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs
 use surf_temp,          only: tsurf_PCM, adapt_tsurf2disappearedice
-use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_tend_h2oice
+use tendencies,         only: compute_tendice, evolve_tend_co2ice, evolve_flux_ssice
 use tracers,            only: adapt_tracers2pressure
 use utility,            only: real2str
@@ -72 +72 @@
 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 :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg/m2]
+real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg/m2]
+real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost  ! Minimum of H2O frost over the last PCM year [kg/m2]
+real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost  ! Minimum of CO2 frost over the last PCM year [kg/m2]
 ! Surface-related:
 real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1 ! Average surface temperature of the second to last PCM run [K]
@@ -141 +141 @@
 
 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,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
+                    q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost,flux_ssice_avg)
 
 ! Initiate soil settings and TI
@@ -169 +169 @@
 call print_msg('> Computing surface ice tendencies',LVL_NFO)
 call compute_tendice(minPCM_h2operice,minPCM_h2ofrost,h2o_ice,d_h2oice)
-call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)),LVL_NFO)
+call print_msg('H2O ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)),LVL_NFO)
 call compute_tendice(minPCM_co2perice,minPCM_co2frost,co2_ice,d_co2ice)
-call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
+call print_msg('CO2 ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
 deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
 
@@ -220 +220 @@
     ! Conversion to surface ice
     call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth)
-    ! We put the sublimating tendency coming from subsurface ice into the overall tendency
-    !where (h2oice_depth > 0. .and. zdqsdif_ssi_tot < -minieps) d_h2oice = zdqsdif_ssi_tot
 end if
 call allocate_loop_state()
@@ -227 +225 @@
 idt = 0
 do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < ntot_yr_sim)
-    call print_msg('**** Iteration of the PEM run [plnt yr]: '//real2str(n_yr_run + dt),LVL_NFO)
+    call print_msg('**** Iteration of the PEM run [plnt y]: '//real2str(n_yr_run + dt),LVL_NFO)
     ! Evolve global surface pressure
     call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg)
@@ -237 +235 @@
     allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
     if (do_layering) then
+        allocate(d_h2oice_new(ngrid,nslope))
         h2oice_depth_old(:,:) = h2oice_depth(:,:)
+        ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency
+        where (h2oice_depth(:,:) > 0. .and. flux_ssice_avg(:,:) < 0._dp) d_h2oice(:,:) = flux_ssice_avg(:,:)
         do islope = 1,nslope
             do i = 1,ngrid
@@ -247 +248 @@
         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)
@@ -253 +253 @@
     else
         zlag(:,:) = 0._dp
+        zshift_surf(:,:) = 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)
@@ -306 +307 @@
             write(num,'("_slope",i2.2)') islope
         end if
-        call write_diagevo('h2oice'//trim(num),'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/))
-        call write_diagevo('co2ice'//trim(num),'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/))
-        call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/))
-        call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/))
+        call write_diagevo('h2oice'//trim(num),'H2O ice','kg/m2',h2o_ice(:,islope),(/dim_ngrid/))
+        call write_diagevo('co2ice'//trim(num),'CO2 ice','kg/m2',co2_ice(:,islope),(/dim_ngrid/))
+        call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg/m2/y',d_h2oice(:,islope),(/dim_ngrid/))
+        call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg/m2/y',d_co2ice(:,islope),(/dim_ngrid/))
         if (co2ice_flow) then
             call write_diagevo('Flow_co2ice'//trim(num),'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
@@ -328 +329 @@
             call write_diagevo('inertiesoil'//trim(num),'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
             if (do_sorption) then
-                call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
-                call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
+                call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg/m2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
+                call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg/m2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
             end if
         end if
@@ -361 +362 @@
     ! Evolve the tendencies
     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_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
-!    else
-!        do i = 1,ngrid
-!            do islope = 1,nslope
-!                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
-!~     end if
+    if (do_layering) then
+        call evolve_flux_ssice(h2oice_depth_old,h2oice_depth,tsurf_avg,tsoil_ts_old,tsoil_ts,flux_ssice_avg)
+    else if (icetable_equilibrium .or. icetable_dynamic) then
+        call evolve_flux_ssice(icetable_depth_old,icetable_depth,tsurf_avg,tsoil_ts_old,tsoil_ts,flux_ssice_avg)
+    end if
 
     ! Increment time
@@ -383 +376 @@
     if (backup_rate > 0) then
         if (mod(idt,backup_rate) == 0) call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
-                                                            icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
+                                                            icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,idt)
     end if
 
@@ -415 +408 @@
 
 call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, &
-                     icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map)
+                     icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map)
 
 ! Update the duration information of the workflow