MODULE tendencies !----------------------------------------------------------------------- ! NAME ! tendencies ! ! DESCRIPTION ! Computation and update of PEM ice evolution tendencies. ! ! AUTHORS & DATE ! R. Vandemeulebrouck ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE compute_tend(ngrid,nslope,min_ice,d_ice) !----------------------------------------------------------------------- ! NAME ! compute_tend ! ! DESCRIPTION ! Compute initial ice evolution tendencies from PCM data. ! ! AUTHORS & DATE ! R. Vandemeulebrouck ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! Based on minima of ice at each point for the PCM years. !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: ngrid integer, intent(in) :: nslope real, dimension(ngrid,nslope,2), intent(in) :: min_ice ! Minima of ice at each point for the PCM years real, dimension(ngrid,nslope), intent(out) :: d_ice ! Evolution of perennial ice ! CODE ! ---- ! We compute the difference d_ice = min_ice(:,:,2) - min_ice(:,:,1) ! If the difference is too small, then there is no evolution where (abs(d_ice) < 1.e-10) d_ice = 0. ! If the minimum over the last year is 0, then we have no perennial ice where (abs(min_ice(:,:,2)) < 1.e-10) d_ice = 0. END SUBROUTINE compute_tend !======================================================================= !======================================================================= SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, & vmr_co2_PCM,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global) !----------------------------------------------------------------------- ! NAME ! recomp_tend_co2 ! ! DESCRIPTION ! Recompute CO2 ice tendency based on pressure and atmospheric changes. ! ! AUTHORS & DATE ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! Adjusts CO2 ice evolution based on Clausius-Clapeyron changes. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: timelen, ngrid, nslope ! Time length, # of grid points and slopes real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM ! CO2 VMR in PCM first layer real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM ! CO2 VMR in PEM first layer real, dimension(ngrid,timelen), intent(in) :: ps_PCM ! Surface pressure in PCM real, intent(in) :: ps_avg_global_ini ! Global average pressure (initial) real, intent(in) :: ps_avg_global ! Global average pressure (current) real, dimension(ngrid,nslope), intent(in) :: d_co2ice_ini ! Initial CO2 ice evolution real, dimension(ngrid,nslope), intent(in) :: co2ice ! CO2 ice surface [kg/m^2] real, dimension(ngrid,nslope), intent(in) :: emissivity ! Emissivity real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! Updated CO2 ice evolution ! LOCAL VARIABLES ! --------------- integer :: i, t, islope real :: coef, avg ! CODE ! ---- write(*,*) "> Updating the CO2 ice tendency for the new pressure" ! Evolution of the water ice for each physical point do i = 1,ngrid do islope = 1,nslope coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2 avg = 0. if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then do t = 1,timelen avg = avg + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM(i,t)/100.)))**4 & - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM(i,t)*(ps_avg_global/ps_avg_global_ini)/100.)))**4 enddo if (avg < 1.e-4) avg = 0. d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*avg/timelen endif enddo enddo END SUBROUTINE recomp_tend_co2 !======================================================================= !======================================================================= SUBROUTINE recomp_tend_h2o(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_PEM_timeseries_old,tsoil_PEM_timeseries_new,d_h2oice) !----------------------------------------------------------------------- ! NAME ! recomp_tend_h2o ! ! DESCRIPTION ! Recompute H2O ice tendency based on soil depth and temperature changes. ! ! AUTHORS & DATE ! JB Clement, 2025 (following E. Vos's work) ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use soil_temp, only: itp_tsoil use subsurface_ice, only: psv ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real, intent(in) :: h2oice_depth_old ! Old H2O ice depth real, intent(in) :: h2oice_depth_new ! New H2O ice depth real, intent(in) :: tsurf ! Surface temperature real, dimension(:,:), intent(in) :: tsoil_PEM_timeseries_old ! Old soil temperature time series real, dimension(:,:), intent(in) :: tsoil_PEM_timeseries_new ! New soil temperature time series real, intent(inout) :: d_h2oice ! Evolution of perennial ice ! LOCAL VARIABLES ! --------------- real :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new integer :: t real, parameter :: coef_diff = 4.e-4 ! Diffusion coefficient real, parameter :: zcdv = 0.0325 ! Drag coefficient ! CODE ! ---- ! Higher resistance due to growing lag layer (higher depth) Rz_old = h2oice_depth_old*zcdv/coef_diff ! Old resistance from PCM Rz_new = h2oice_depth_new*zcdv/coef_diff ! New resistance based on new depth R_dec = Rz_old/Rz_new ! Decrease because of resistance ! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location psv_max_old = 0. psv_max_new = 0. do t = 1,size(tsoil_PEM_timeseries_old,2) psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_PEM_timeseries_old(:,t),tsurf,h2oice_depth_old))) psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_PEM_timeseries_new(:,t),tsurf,h2oice_depth_new))) enddo ! Lower humidity due to growing lag layer (higher depth) if (abs(psv_max_old) < 1.e2*epsilon(1.)) then hum_dec = 1. else hum_dec = psv_max_new/psv_max_old ! Decrease because of lower water vapor pressure at the new depth endif ! Flux correction (decrease) d_h2oice = d_h2oice*R_dec*hum_dec END SUBROUTINE recomp_tend_h2o !======================================================================= END MODULE tendencies