module co2glaciers_mod implicit none LOGICAL co2glaciersflow ! 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_slope,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_slope(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) ! Physical field: CO2 condensation temperature [K] REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2 glacier before flow !----------------------------- call computeTcond(timelen,ngrid,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) call compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax) call transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh) RETURN end subroutine subroutine compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax) USE comconst_mod, ONLY: pi,g !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Compute the maximum thickness of CO2 glaciers given a slope angle !!! before initating flow !!! !!! Author: LL, based on theoretical work by A.Grau Galofre (LPG) !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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) :: Tcond(ngrid) ! Physical field: CO2 condensation temperature [K] REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg] ! Outputs REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum co2 thickness before flaw [m] ! Local INTEGER,PARAMETER :: n = 7 ! flow law exponent Nye et al., 2000 REAL,PARAMETER :: Rg = 8.3145 ! gas constant [J/K/mol] REAL,PARAMETER :: Q = 59000. ! Activation Energy [J/mol], Nye et al., 2000 DOUBLE PRECISION,PARAMETER :: C = 1.8138e-21 ! Nye et al., 2000 [s/m^n] DOUBLE PRECISION :: Ad = 1.202e11 ! Softness prefactor [MPa^-n] Nye et al., 2000 REAL :: Ro,Ho, S,ratio ! gemoetry from Nye et al., 2000 DOUBLE PRECISION :: A,Ao ! softness parameter [s/m^n] DOUBLE PRECISION :: C1 ! intermediate variable DOUBLE PRECISION :: t_0 ! relaxation time (assuming radial symetry) [s] DOUBLE PRECISION :: u ! characteristic horizontal deformation rate [m/s] DOUBLE PRECISION :: tau_d ! characteristic basal drag, understood as the strest that an ice CO2 mass flowing under its weight balanced by viscosity REAL :: rho_co2(ngrid) ! co2 ice density [kg/m^3] INTEGER :: ig,islope ! loop variables REAL :: slo_angle ! 0. Geometry parameters Ro = 200e3 Ho = 3000. ratio = 2./3. S = Ho/Ro*1/((2+1./n)*(1+1./n))*(1-ratio**(1+1./n))**(1./(2+1./n)-1)*ratio**(1+1./n-1) ! 1. Flow parameters Ao = 3**(1./(2*n+2))*Ad do ig = 1,n Ao = Ao*1e-6 enddo ! 2. Compute rho_co2 DO ig = 1,ngrid rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3 ! Mangan et al. 2017 ENDDO ! 3. Compute max thickness DO ig = 1,ngrid A = Ao*exp(-Q/(Rg*Tcond(ig))) C1 = A*(rho_co2(ig)*g)**(float(n))/float(n+2) t_0 = Ro/(C1*(5*n+3))*(float(2*n+1)/float(n+1))**n u = Ro/t_0 tau_d = (u*Ho*(rho_co2(ig)*g)**2*S**2/(2*A))**(1./(n+2)) DO islope = 1,nslope if(islope.eq.iflat) then hmax(ig,islope) = 1.e6 else slo_angle = abs(def_slope_mean(islope)*pi/180.) hmax(ig,islope) = tau_d/(rho_co2(ig)*g*slo_angle) endif ENDDO ENDDO RETURN end subroutine subroutine transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! !!! Purpose: Transfer the excess of ice from one subslope to another !!! No transfer between mesh at the time !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! USE comconst_mod, ONLY: pi 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 CO2 glaciers before initiating flow [m] REAL, INTENT(IN) :: Tcond(ngrid) ! CO2 condensation temperature [K] ! Outputs REAL, INTENT(INOUT) :: co2ice_slope(ngrid,nslope) ! CO2 in the subslope [kg/m^2] REAL, INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope REAL, INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! boolean to check if there is flow in the mesh ! Local INTEGER ig,islope ! loop REAL rho_co2(ngrid) ! density of CO2, temperature dependant [kg/m^3] INTEGER iaval ! ice will be transfered here ! 0. Compute rho_co2 DO ig = 1,ngrid rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3 ! Mangan et al. 2017 ENDDO ! 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(co2ice_slope(ig,islope).ge.rho_co2(ig)*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 co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + & (co2ice_slope(ig,islope) - rho_co2(ig)*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.) co2ice_slope(ig,islope)=rho_co2(ig)*hmax(ig,islope) * & cos(pi*def_slope_mean(islope)/180.) flag_co2flow(ig,islope) = 1. flag_co2flow_mesh(ig) = 1. ENDIF ! co2ice > hmax ENDIF ! iflat ENDDO !islope ENDDO !ig RETURN end subroutine subroutine computeTcond(timelen,ngrid,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 ! # 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) ! 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_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100)) ENDDO Tcond(ig) = ave/timelen ENDDO RETURN end subroutine end module