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 ! 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 ! 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 glacier" 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,intent(in) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles real,intent(in) :: Tice(ngrid,nslope) ! Ice Temperature [K] ! Ouputs: real,intent(inout) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] real,intent(inout) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes real,intent(inout) :: flag_h2oflow_mesh(ngrid) ! same but within the mesh ! Local real :: hmax(ngrid,nslope) ! 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 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 ! arguments ! -------- ! Inputs integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes integer,intent(in) :: iflat ! index of the flat subslope real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg] real,intent(in) :: Tice(ngrid,nslope) ! Physical field: ice temperature [K] character(len=3), intent(in) :: name_ice ! Nature of the ice ! Outputs real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum thickness before flaw [m] ! Local DOUBLE PRECISION :: 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 real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3] integer :: ig,islope ! loop variables real :: slo_angle ! 1. Compute rho if(name_ice.eq."co2") then DO ig = 1,ngrid DO islope = 1,nslope rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3 ! Mangan et al. 2017 tau_d = 5.e3 ENDDO ENDDO elseif (name_ice.eq."h2o") then DO ig = 1,ngrid DO islope = 1,nslope rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 tau_d = 1.e5 ENDDO ENDDO else call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1) endif ! 3. Compute max thickness DO ig = 1,ngrid DO islope = 1,nslope if(islope.eq.iflat) then hmax(ig,islope) = 1.e8 else slo_angle = abs(def_slope_mean(islope)*pi/180.) hmax(ig,islope) = tau_d/(rho(ig,islope)*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 abort_pem_mod, only: abort_pem #ifndef CPP_STD use comcstfi_h, only: pi #else use comcstfi_mod, only: pi #endif implicit none ! arguments ! -------- ! Inputs integer, intent(in) :: ngrid,nslope !# of physical points and subslope integer, intent(in) :: iflat ! index of the flat subslope real, intent(in) :: subslope_dist(ngrid,nslope) ! Distribution of the subgrid slopes within the mesh real, intent(in) :: def_slope_mean(nslope) ! values of the subgrid slopes real, intent(in) :: hmax(ngrid,nslope) ! maximum height of the glaciers before initiating flow [m] real, intent(in) :: Tice(ngrid,nslope) ! Ice temperature[K] character(len=3), intent(in) :: name_ice ! Nature of the ice ! Outputs real, intent(inout) :: qice(ngrid,nslope) ! CO2 in the subslope [kg/m^2] real, intent(inout) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope real, intent(inout) :: flag_flowmesh(ngrid) ! boolean to check if there is flow in the mesh ! Local integer ig,islope ! loop real rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3] integer iaval ! ice will be transfered here ! 0. Compute rho if(name_ice.eq."co2") then DO ig = 1,ngrid DO islope = 1,nslope rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3 ! Mangan et al. 2017 ENDDO ENDDO elseif (name_ice.eq."h2o") then DO ig = 1,ngrid DO islope = 1,nslope rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 ENDDO ENDDO else call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1) endif ! 1. Compute the transfer of ice DO ig = 1,ngrid DO islope = 1,nslope IF(islope.ne.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).ge.rho(ig,islope)*hmax(ig,islope) * & cos(pi*def_slope_mean(islope)/180.)) THEN ! Second: determine the flatest slopes possible: IF(islope.gt.iflat) THEN iaval=islope-1 ELSE iaval=islope+1 ENDIF do while ((iaval.ne.iflat).and. & (subslope_dist(ig,iaval).eq.0)) IF(iaval.gt.iflat) THEN iaval=iaval-1 ELSE iaval=iaval+1 ENDIF enddo qice(ig,iaval) = qice(ig,iaval) + & (qice(ig,islope) - rho(ig,islope)*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(ig,islope)*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