MODULE tendencies !----------------------------------------------------------------------- ! NAME ! tendencies ! ! DESCRIPTION ! Computation and update of PEM ice evolution tendencies. ! ! AUTHORS & DATE ! R. Vandemeulebrouck ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use numerics, only: dp, di, k4, minieps, tol ! DECLARATION ! ----------- implicit none contains !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE compute_tendice(min_ice,is_perice,d_ice) !----------------------------------------------------------------------- ! NAME ! compute_tendice ! ! 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 ! --------- real(dp), dimension(:,:,:), intent(in) :: min_ice logical(k4), dimension(:,:), intent(in) :: is_perice real(dp), dimension(:,:), intent(out) :: d_ice ! CODE ! ---- ! We compute the difference to get the tendency d_ice = min_ice(:,:,2) - min_ice(:,:,1) ! If the difference is too small, then there is no evolution where (abs(d_ice) < minieps) d_ice = 0._dp ! If the tendency is negative but there is no ice reservoir for the PEM where (abs(d_ice) < 0._dp .and. .not. is_perice) d_ice = 0._dp END SUBROUTINE compute_tendice !======================================================================= !======================================================================= SUBROUTINE evolve_tend_co2ice(d_co2ice_ini,co2ice,emissivity,q_co2_ts_ini,q_co2_ts,ps_PCM,ps_avg_glob_ini,ps_avg_global,d_co2ice) !----------------------------------------------------------------------- ! NAME ! evolve_tend_co2ice ! ! 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 geometry, only: ngrid, nslope, nday use planet, only: alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, m_co2, A, B use orbit, only: yr_len_sol, sol_len_s use display, only: print_msg use utility, only: real2str ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), dimension(:,:), intent(in) :: q_co2_ts, q_co2_ts_ini, ps_PCM, d_co2ice_ini, co2ice, emissivity real(dp), intent(in) :: ps_avg_global, ps_avg_glob_ini real(dp), dimension(:,:), intent(inout) :: d_co2ice ! LOCAL VARIABLES ! --------------- integer(di) :: i, iday, islope real(dp) :: coef, avg real(dp), dimension(ngrid,nday) :: vmr_co2_ts, vmr_co2_ts_ini ! CODE ! ---- call print_msg("> Evolving the CO2 ice tendency according to the new pressure") call print_msg('Old CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice))) ! Compute volume mixing ratio vmr_co2_ts(:,:) = q_co2_ts(:,:)/(A*q_co2_ts(:,:) + B)/m_co2 vmr_co2_ts_ini(:,:) = q_co2_ts_ini(:,:)/(A*q_co2_ts_ini(:,:) + B)/m_co2 ! Recompute the tendency do i = 1,ngrid do islope = 1,nslope coef = yr_len_sol*sol_len_s*emissivity(i,islope)*sigmaB/Lco2 avg = 0._dp if (co2ice(i,islope) > 1.e-4_dp .and. abs(d_co2ice(i,islope)) > 1.e-5_dp) then do iday = 1,nday avg = avg + (beta_clap_co2/(alpha_clap_co2 - log(q_co2_ts_ini(i,iday)*ps_PCM(i,iday)/100._dp)))**4 & - (beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_ts(i,iday)*ps_PCM(i,iday)*(ps_avg_global/ps_avg_glob_ini)/100._dp)))**4 end do if (avg < 1.e-4_dp) avg = 0._dp d_co2ice(i,islope) = d_co2ice_ini(i,islope) - coef*avg/nday end if end do end do call print_msg('New CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str( minval(d_co2ice))//' | '//real2str(maxval(d_co2ice))) END SUBROUTINE evolve_tend_co2ice !======================================================================= !======================================================================= SUBROUTINE evolve_tend_h2oice(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_ts_old,tsoil_ts_new,d_h2oice) !----------------------------------------------------------------------- ! NAME ! evolve_tend_h2oice ! ! 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 use display, only: print_msg ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- real(dp), intent(in) :: h2oice_depth_old ! Old H2O ice depth real(dp), intent(in) :: h2oice_depth_new ! New H2O ice depth real(dp), intent(in) :: tsurf ! Surface temperature real(dp), dimension(:,:), intent(in) :: tsoil_ts_old ! Old soil temperature time series real(dp), dimension(:,:), intent(in) :: tsoil_ts_new ! New soil temperature time series real(dp), intent(inout) :: d_h2oice ! Evolution of perennial ice ! LOCAL VARIABLES ! --------------- real(dp) :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new integer(di) :: iday real(dp), parameter :: coef_diff = 4.e-4 ! Diffusion coefficient real(dp), parameter :: zcdv = 0.0325 ! Drag coefficient ! CODE ! ---- call print_msg("> Updating the H2O sub-surface ice tendency due to lag layer") ! 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._dp psv_max_new = 0._dp do iday = 1,size(tsoil_ts_old,2) psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_ts_old(:,iday),tsurf,h2oice_depth_old))) psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_ts_new(:,iday),tsurf,h2oice_depth_new))) end do ! Lower humidity due to growing lag layer (higher depth) if (abs(psv_max_old) < tol) then hum_dec = 1._dp else hum_dec = psv_max_new/psv_max_old ! Decrease because of lower water vapor pressure at the new depth end if ! Flux correction (decrease) d_h2oice = d_h2oice*R_dec*hum_dec END SUBROUTINE evolve_tend_h2oice !======================================================================= END MODULE tendencies