module glaciers_mod implicit none LOGICAL co2glaciersflow ! True by default, to compute co2 ice flow. Read in pem.def LOGICAL h2oglaciersflow ! True by default, to compute co2 ice flow. Read in pem.def !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Compute CO2 glacier flows !!! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! contains subroutine co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_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,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) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol] REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical x Time field: surface pressure given by the GCM [Pa] REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surface pressure from the GCM [Pa] REAL,INTENT(IN) :: global_ave_ps_PEM ! global averaged surface pressure during the PEM iteration [Pa] ! Ouputs: REAL,INTENT(INOUT) :: co2ice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! same but within the mesh ! Local REAL :: Tcond(ngrid,nslope) ! Physical field: CO2 condensation temperature [K] REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2 glacier before flow !----------------------------- call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_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) RETURN end subroutine subroutine h2oglaciers_evol(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 !----------------------------- 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) RETURN end subroutine 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 RETURN end subroutine 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 RETURN end subroutine subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_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,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol] REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa] REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa] REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration ! OUTPUT REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged ! LOCAL INTEGER :: ig,it,islope ! 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_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100)) ENDDO DO islope = 1,nslope Tcond(ig,islope) = ave/timelen ENDDO ENDDO RETURN end subroutine end module