MODULE glaciers_mod implicit none ! Flags for ice management logical :: h2oice_flow ! True by default, to compute H2O ice flow. Read in "run_PEM.def" logical :: co2ice_flow ! True by default, to compute CO2 ice flow. Read in "run_PEM.def" logical :: metam_h2oice ! False by default, to compute H2O ice metamorphism. Read in "run_PEM.def" logical :: metam_co2ice ! False by default, to compute CO2 ice metamorphism. Read in "run_PEM.def" ! Thresholds for ice management real, save :: inf_h2oice_threshold ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2] real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice [kg.m-2] real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice [kg.m-2] !======================================================================= contains !======================================================================= SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do !!! the ice transfer !!! !!! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none ! Inputs: !-------- 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 ! Physical points x Slopes: Distribution of the subgrid slopes real, dimension(ngrid), intent(in) :: def_slope_mean ! Physical points: values of the sub grid slope angles real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM ! Physical x Time field : VMR of co2 in the first layer [mol/mol] real, dimension(ngrid,timelen), intent(in) :: ps_PCM ! Physical x Time field: surface pressure given by the PCM [Pa] real, intent(in) :: global_avg_ps_PCM ! Global averaged surface pressure from the PCM [Pa] real, intent(in) :: global_avg_ps_PEM ! global averaged surface pressure during the PEM iteration [Pa] ! Ouputs: !-------- real, dimension(ngrid,nslope), intent(inout) :: co2ice ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow ! flag to see if there is flow on the subgrid slopes real, dimension(ngrid), intent(inout) :: flag_co2flow_mesh ! same but within the mesh ! Local !------ real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K] real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2 glacier before flow write(*,*) "Flow of CO2 glaciers" call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,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,"co2",co2ice,flag_co2flow,flag_co2flow_mesh) END SUBROUTINE flow_co2glaciers !======================================================================= SUBROUTINE flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do !!! the ice transfer !!! !!! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none ! arguments ! --------- ! Inputs: 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 ! Physical points x Slopes : Distribution of the subgrid slopes real, dimension(ngrid), intent(in) :: def_slope_mean ! Slopes: values of the sub grid slope angles real, dimension(ngrid,nslope), intent(in) :: Tice ! Ice Temperature [K] ! Ouputs: real, dimension(ngrid,nslope), intent(inout) :: h2oice ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] real, dimension(ngrid,nslope), intent(inout) :: flag_h2oflow ! flag to see if there is flow on the subgrid slopes real, dimension(ngrid), intent(inout) :: flag_h2oflow_mesh ! same but within the mesh ! Local real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2 glacier before flow 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,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh) END SUBROUTINE flow_h2oglaciers !======================================================================= SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax) use ice_table_mod, only: rho_ice use abort_pem_mod, only: abort_pem #ifndef CPP_STD use comcstfi_h, only: pi, g #else use comcstfi_mod, only: pi, g #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow !!! !!! Author: LL,based on work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022) !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none ! Inputs ! ------ 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 ! Slope field: Values of the subgrid slope angles [deg] real, dimension(ngrid,nslope), intent(in) :: Tice ! Physical field: ice temperature [K] character(3), intent(in) :: name_ice ! Nature of ice ! Outputs ! ------- real, dimension(ngrid,nslope), intent(out) :: hmax ! Physical grid x Slope field: maximum thickness before flaw [m] ! Local ! ----- 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 ! loop variables real :: slo_angle select case (trim(adjustl(name_ice))) case('h2o') tau_d = 1.e5 case('co2') tau_d = 5.e3 case default call abort_pem("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,name_ice,qice,flag_flow,flag_flowmesh) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Transfer the excess of ice from one subslope to another !!! No transfer between mesh at the time !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use ice_table_mod, only: rho_ice use abort_pem_mod, only: abort_pem #ifndef CPP_STD use comcstfi_h, only: pi #else use comcstfi_mod, only: pi #endif implicit none ! Inputs !------- 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 within the mesh real, dimension(nslope), intent(in) :: def_slope_mean ! values of the subgrid slopes real, dimension(ngrid,nslope), intent(in) :: hmax ! maximum height of the glaciers before initiating flow [m] real, dimension(ngrid,nslope), intent(in) :: Tice ! Ice temperature[K] character(3), intent(in) :: name_ice ! Nature of the ice ! Outputs !-------- real, dimension(ngrid,nslope), intent(inout) :: qice ! CO2 in the subslope [kg/m^2] real, dimension(ngrid,nslope), intent(inout) :: flag_flow ! boolean to check if there is flow on a subgrid slope real, dimension(ngrid), intent(inout) :: flag_flowmesh ! boolean to check if there is flow in the mesh ! Local !------ integer :: ig, islope ! loop integer :: iaval ! ice will be transfered here 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. flag_flowmesh(ig) = 1. endif ! co2ice > hmax endif ! iflat enddo !islope enddo !ig END SUBROUTINE !======================================================================= SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Compute CO2 condensation temperature !!! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2 implicit none ! arguments: ! ---------- ! INPUT integer, intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM ! Physical points x times field: VMR of CO2 in the first layer [mol/mol] real, dimension(ngrid,timelen), intent(in) :: ps_PCM ! Physical points x times field: surface pressure in the PCM [Pa] real, intent(in) :: global_avg_ps_PCM ! Global averaged surfacepressure in the PCM [Pa] real, intent(in) :: global_avg_ps_PEM ! Global averaged surface pressure computed during the PEM iteration ! OUTPUT real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged ! LOCAL integer :: ig, it ! For loop real :: ave ! Intermediate to compute average do ig = 1,ngrid ave = 0 do it = 1,timelen ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*global_avg_ps_PCM/global_avg_ps_PEM/100)) enddo Tcond(ig,:) = ave/timelen enddo END SUBROUTINE computeTcondCO2 END MODULE glaciers_mod