MODULE glaciers !----------------------------------------------------------------------- ! NAME ! glaciers ! ! DESCRIPTION ! Compute flow and transfer of CO2 and H2O ice glaciers on slopes ! based on maximum ice thickness and basal drag conditions. ! ! AUTHORS & DATE ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! PARAMETERS ! ---------- logical :: h2oice_flow ! True by default, flag to compute H2O ice flow logical :: co2ice_flow ! True by default, flag to compute CO2 ice flow contains !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !======================================================================= SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,co2ice,flag_co2flow) !----------------------------------------------------------------------- ! NAME ! flow_co2glaciers ! ! DESCRIPTION ! Compute maximum thickness and transfer CO2 ice between subslopes. ! ! AUTHORS & DATE ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope real, dimension(ngrid,nslope), intent(in) :: subslope_dist ! Grid points x Slopes: Distribution of the subgrid slopes real, dimension(ngrid), intent(in) :: def_slope_mean ! Grid points: values of the sub grid slope angles real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM ! Grid points x Time field : VMR of CO2 in the first layer [mol/mol] real, dimension(ngrid,timelen), intent(in) :: ps_PCM ! Grid points x Time field: surface pressure given by the PCM [Pa] real, intent(in) :: ps_avg_global_ini ! Global averaged surface pressure at the beginning [Pa] real, intent(in) :: ps_avg_global ! Global averaged surface pressure during the PEM iteration [Pa] real, dimension(ngrid,nslope), intent(inout) :: co2ice ! Grid points x Slope field: CO2 ice on the subgrid slopes [kg/m^2] integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_co2flow ! Flag to see if there is flow on the subgrid slopes ! LOCAL VARIABLES ! --------------- real, dimension(ngrid,nslope) :: Tcond ! CO2 condensation temperature [K] real, dimension(ngrid,nslope) :: hmax ! Maximum thickness before flow ! CODE ! ---- write(*,*) "> Flow of CO2 glaciers" call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond) call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax) call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tcond,co2ice,flag_co2flow) END SUBROUTINE flow_co2glaciers !======================================================================= !======================================================================= SUBROUTINE flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow) !----------------------------------------------------------------------- ! NAME ! flow_h2oglaciers ! ! DESCRIPTION ! Compute maximum thickness and transfer H2O ice between subslopes. ! ! AUTHORS & DATE ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope real, dimension(ngrid,nslope), intent(in) :: subslope_dist ! Distribution of the subgrid slopes real, dimension(ngrid), intent(in) :: def_slope_mean ! Values of the sub grid slope angles real, dimension(ngrid,nslope), intent(in) :: Tice ! Ice Temperature [K] real, dimension(ngrid,nslope), intent(inout) :: h2oice ! H2O ice on the subgrid slopes [kg/m^2] integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_h2oflow ! Flow flag on subgrid slopes ! LOCAL VARIABLES ! --------------- real, dimension(ngrid,nslope) :: hmax ! Maximum thickness before flow ! CODE ! ---- write(*,*) "> Flow of H2O glaciers" call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax) call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,h2oice,flag_h2oflow) END SUBROUTINE flow_h2oglaciers !======================================================================= !======================================================================= SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax) !----------------------------------------------------------------------- ! NAME ! compute_hmaxglaciers ! ! DESCRIPTION ! Compute the maximum thickness of CO2 and H2O glaciers given a ! slope angle before initiating flow. ! ! AUTHORS & DATE ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! Based on work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022) ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use ice_table, only: rho_ice use stop_pem, only: stop_clean use phys_constants, only: pi, g ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: ngrid, nslope ! # of grid points and subslopes integer, intent(in) :: iflat ! index of the flat subslope real, dimension(nslope), intent(in) :: def_slope_mean ! Values of the subgrid slope angles [deg] real, dimension(ngrid,nslope), intent(in) :: Tice ! Ice temperature [K] character(3), intent(in) :: name_ice ! Nature of ice real, dimension(ngrid,nslope), intent(out) :: hmax ! Maximum thickness before flow [m] ! LOCAL VARIABLES ! --------------- real :: tau_d ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith integer :: ig, islope real :: slo_angle ! CODE ! ---- select case (trim(adjustl(name_ice))) case('h2o') tau_d = 1.e5 case('co2') tau_d = 5.e3 case default call stop_clean("compute_hmaxglaciers","Type of ice not known!",1) end select do ig = 1,ngrid do islope = 1,nslope if (islope == iflat) then hmax(ig,islope) = 1.e8 else slo_angle = abs(def_slope_mean(islope)*pi/180.) hmax(ig,islope) = tau_d/(rho_ice(Tice(ig,islope),name_ice)*g*slo_angle) endif enddo enddo END SUBROUTINE compute_hmaxglaciers !======================================================================= !======================================================================= SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,qice,flag_flow) !----------------------------------------------------------------------- ! NAME ! transfer_ice_duringflow ! ! DESCRIPTION ! Transfer the excess of ice from one subslope to another. ! No transfer between mesh at the time. ! ! AUTHORS & DATE ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use ice_table, only: rho_ice use stop_pem, only: stop_clean use phys_constants, only: pi ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: ngrid, nslope ! # of physical points and subslope integer, intent(in) :: iflat ! Index of the flat subslope real, dimension(ngrid,nslope), intent(in) :: subslope_dist ! Distribution of the subgrid slopes real, dimension(nslope), intent(in) :: def_slope_mean ! Values of the subgrid slopes real, dimension(ngrid,nslope), intent(in) :: hmax ! Maximum height before initiating flow [m] real, dimension(ngrid,nslope), intent(in) :: Tice ! Ice temperature [K] real, dimension(ngrid,nslope), intent(inout) :: qice ! Ice in the subslope [kg/m^2] integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_flow ! Flow flag on subgrid slopes ! LOCAL VARIABLES ! --------------- integer :: ig, islope integer :: iaval ! index where ice is transferred ! CODE ! ---- flag_flow = 0 do ig = 1,ngrid do islope = 1,nslope if (islope /= iflat) then ! ice can be infinite on flat ground ! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely on flat ground if (qice(ig,islope) >= rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) then ! Second: determine the flatest slopes possible: if (islope > iflat) then iaval = islope - 1 else iaval = islope + 1 endif do while (iaval /= iflat .and. subslope_dist(ig,iaval) == 0) if (iaval > iflat) then iaval = iaval - 1 else iaval = iaval + 1 endif enddo qice(ig,iaval) = qice(ig,iaval) + (qice(ig,islope) - rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.)) & *subslope_dist(ig,islope)/subslope_dist(ig,iaval)*cos(pi*def_slope_mean(iaval)/180.)/cos(pi*def_slope_mean(islope)/180.) qice(ig,islope) = rho_ice(Tice(ig,islope),'h2o')*hmax(ig,islope)*cos(pi*def_slope_mean(islope)/180.) flag_flow(ig,islope) = 1 endif ! co2ice > hmax endif ! iflat enddo !islope enddo !ig END SUBROUTINE !======================================================================= SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond) !----------------------------------------------------------------------- ! NAME ! computeTcondCO2 ! ! DESCRIPTION ! Compute CO2 condensation temperature. ! ! AUTHORS & DATE ! L. Lange ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2 ! DECLARATION ! ----------- implicit none ! ARGUMENTS ! --------- integer, intent(in) :: timelen, ngrid, nslope real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM ! VMR of CO2 in the first layer [mol/mol] real, dimension(ngrid,timelen), intent(in) :: ps_PCM ! Surface pressure in the PCM [Pa] real, intent(in) :: ps_avg_global_ini ! Global averaged surfacepressure in the PCM [Pa] real, intent(in) :: ps_avg_global ! Global averaged surface pressure computed during the PEM iteration real, dimension(ngrid,nslope), intent(out) :: Tcond ! Condensation temperature of CO2, yearly averaged ! LOCAL VARIABLES ! --------------- integer :: ig, it ! CODE ! ---- do ig = 1,ngrid Tcond(ig,:) = 0 do it = 1,timelen Tcond(ig,:) = Tcond(ig,:) + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*ps_avg_global_ini/ps_avg_global/100)) enddo Tcond(ig,:) = Tcond(ig,:)/timelen enddo END SUBROUTINE computeTcondCO2 !======================================================================= END MODULE glaciers