| 1 | module co2glaciers_mod |
|---|
| 2 | implicit none |
|---|
| 3 | LOGICAL co2glaciersflow ! True by default, to compute co2 ice flow. Read in pem.def |
|---|
| 4 | |
|---|
| 5 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 6 | !!! |
|---|
| 7 | !!! Purpose: Compute CO2 glacier flows |
|---|
| 8 | !!! |
|---|
| 9 | !!! Author: LL |
|---|
| 10 | !!! |
|---|
| 11 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 12 | |
|---|
| 13 | contains |
|---|
| 14 | |
|---|
| 15 | |
|---|
| 16 | 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) |
|---|
| 17 | |
|---|
| 18 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 19 | !!! |
|---|
| 20 | !!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do |
|---|
| 21 | !!! the ice transfer |
|---|
| 22 | !!! |
|---|
| 23 | !!! |
|---|
| 24 | !!! Author: LL |
|---|
| 25 | !!! |
|---|
| 26 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 27 | |
|---|
| 28 | |
|---|
| 29 | IMPLICIT NONE |
|---|
| 30 | |
|---|
| 31 | ! arguments |
|---|
| 32 | ! --------- |
|---|
| 33 | |
|---|
| 34 | ! Inputs: |
|---|
| 35 | INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! # number of time sample, physical points, subslopes, index of the flat subslope |
|---|
| 36 | 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 |
|---|
| 37 | REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol] |
|---|
| 38 | REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical x Time field: surface pressure given by the GCM [Pa] |
|---|
| 39 | REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surface pressure from the GCM [Pa] |
|---|
| 40 | REAL,INTENT(IN) :: global_ave_ps_PEM ! global averaged surface pressure during the PEM iteration [Pa] |
|---|
| 41 | |
|---|
| 42 | ! Ouputs: |
|---|
| 43 | REAL,INTENT(INOUT) :: co2ice_slope(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] |
|---|
| 44 | REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes |
|---|
| 45 | REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! same but within the mesh |
|---|
| 46 | |
|---|
| 47 | |
|---|
| 48 | ! Local |
|---|
| 49 | REAL :: Tcond(ngrid) ! Physical field: CO2 condensation temperature [K] |
|---|
| 50 | REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2 glacier before flow |
|---|
| 51 | |
|---|
| 52 | !----------------------------- |
|---|
| 53 | call computeTcond(timelen,ngrid,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) |
|---|
| 54 | |
|---|
| 55 | call compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax) |
|---|
| 56 | |
|---|
| 57 | call transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh) |
|---|
| 58 | RETURN |
|---|
| 59 | end subroutine |
|---|
| 60 | |
|---|
| 61 | |
|---|
| 62 | |
|---|
| 63 | subroutine compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax) |
|---|
| 64 | |
|---|
| 65 | USE comconst_mod, ONLY: pi,g |
|---|
| 66 | |
|---|
| 67 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 68 | !!! |
|---|
| 69 | !!! Purpose: Compute the maximum thickness of CO2 glaciers given a slope angle |
|---|
| 70 | !!! before initating flow |
|---|
| 71 | !!! |
|---|
| 72 | !!! Author: LL, based on theoretical work by A.Grau Galofre (LPG) |
|---|
| 73 | !!! |
|---|
| 74 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 75 | |
|---|
| 76 | IMPLICIT NONE |
|---|
| 77 | |
|---|
| 78 | ! arguments |
|---|
| 79 | ! -------- |
|---|
| 80 | |
|---|
| 81 | ! Inputs |
|---|
| 82 | INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes |
|---|
| 83 | INTEGER,INTENT(IN) :: iflat ! index of the flat subslope |
|---|
| 84 | REAL,INTENT(IN) :: Tcond(ngrid) ! Physical field: CO2 condensation temperature [K] |
|---|
| 85 | REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg] |
|---|
| 86 | ! Outputs |
|---|
| 87 | REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum co2 thickness before flaw [m] |
|---|
| 88 | ! Local |
|---|
| 89 | INTEGER,PARAMETER :: n = 7 ! flow law exponent Nye et al., 2000 |
|---|
| 90 | REAL,PARAMETER :: Rg = 8.3145 ! gas constant [J/K/mol] |
|---|
| 91 | REAL,PARAMETER :: Q = 59000. ! Activation Energy [J/mol], Nye et al., 2000 |
|---|
| 92 | DOUBLE PRECISION,PARAMETER :: C = 1.8138e-21 ! Nye et al., 2000 [s/m^n] |
|---|
| 93 | DOUBLE PRECISION :: Ad = 1.202e11 ! Softness prefactor [MPa^-n] Nye et al., 2000 |
|---|
| 94 | REAL :: Ro,Ho, S,ratio ! gemoetry from Nye et al., 2000 |
|---|
| 95 | DOUBLE PRECISION :: A,Ao ! softness parameter [s/m^n] |
|---|
| 96 | DOUBLE PRECISION :: C1 ! intermediate variable |
|---|
| 97 | DOUBLE PRECISION :: t_0 ! relaxation time (assuming radial symetry) [s] |
|---|
| 98 | DOUBLE PRECISION :: u ! characteristic horizontal deformation rate [m/s] |
|---|
| 99 | DOUBLE PRECISION :: tau_d ! characteristic basal drag, understood as the strest that an ice CO2 mass flowing under its weight balanced by viscosity |
|---|
| 100 | REAL :: rho_co2(ngrid) ! co2 ice density [kg/m^3] |
|---|
| 101 | INTEGER :: ig,islope ! loop variables |
|---|
| 102 | REAL :: slo_angle |
|---|
| 103 | |
|---|
| 104 | ! 0. Geometry parameters |
|---|
| 105 | Ro = 200e3 |
|---|
| 106 | Ho = 3000. |
|---|
| 107 | ratio = 2./3. |
|---|
| 108 | 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) |
|---|
| 109 | ! 1. Flow parameters |
|---|
| 110 | Ao = 3**(1./(2*n+2))*Ad |
|---|
| 111 | do ig = 1,n |
|---|
| 112 | Ao = Ao*1e-6 |
|---|
| 113 | enddo |
|---|
| 114 | ! 2. Compute rho_co2 |
|---|
| 115 | DO ig = 1,ngrid |
|---|
| 116 | rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3 ! Mangan et al. 2017 |
|---|
| 117 | ENDDO |
|---|
| 118 | ! 3. Compute max thickness |
|---|
| 119 | DO ig = 1,ngrid |
|---|
| 120 | A = Ao*exp(-Q/(Rg*Tcond(ig))) |
|---|
| 121 | C1 = A*(rho_co2(ig)*g)**(float(n))/float(n+2) |
|---|
| 122 | t_0 = Ro/(C1*(5*n+3))*(float(2*n+1)/float(n+1))**n |
|---|
| 123 | u = Ro/t_0 |
|---|
| 124 | tau_d = (u*Ho*(rho_co2(ig)*g)**2*S**2/(2*A))**(1./(n+2)) |
|---|
| 125 | DO islope = 1,nslope |
|---|
| 126 | if(islope.eq.iflat) then |
|---|
| 127 | hmax(ig,islope) = 1.e6 |
|---|
| 128 | else |
|---|
| 129 | slo_angle = abs(def_slope_mean(islope)*pi/180.) |
|---|
| 130 | hmax(ig,islope) = tau_d/(rho_co2(ig)*g*slo_angle) |
|---|
| 131 | endif |
|---|
| 132 | ENDDO |
|---|
| 133 | ENDDO |
|---|
| 134 | RETURN |
|---|
| 135 | |
|---|
| 136 | end subroutine |
|---|
| 137 | |
|---|
| 138 | subroutine transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh) |
|---|
| 139 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 140 | !!! |
|---|
| 141 | !!! Purpose: Transfer the excess of ice from one subslope to another |
|---|
| 142 | !!! No transfer between mesh at the time |
|---|
| 143 | !!! Author: LL |
|---|
| 144 | !!! |
|---|
| 145 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 146 | |
|---|
| 147 | USE comconst_mod, ONLY: pi |
|---|
| 148 | |
|---|
| 149 | |
|---|
| 150 | implicit none |
|---|
| 151 | |
|---|
| 152 | ! arguments |
|---|
| 153 | ! -------- |
|---|
| 154 | |
|---|
| 155 | ! Inputs |
|---|
| 156 | INTEGER, INTENT(IN) :: ngrid,nslope !# of physical points and subslope |
|---|
| 157 | INTEGER, INTENT(IN) :: iflat ! index of the flat subslope |
|---|
| 158 | REAL, INTENT(IN) :: subslope_dist(ngrid,nslope) ! Distribution of the subgrid slopes within the mesh |
|---|
| 159 | REAL, INTENT(IN) :: def_slope_mean(nslope) ! values of the subgrid slopes |
|---|
| 160 | REAL, INTENT(IN) :: hmax(ngrid,nslope) ! maximum height of the CO2 glaciers before initiating flow [m] |
|---|
| 161 | REAL, INTENT(IN) :: Tcond(ngrid) ! CO2 condensation temperature [K] |
|---|
| 162 | ! Outputs |
|---|
| 163 | REAL, INTENT(INOUT) :: co2ice_slope(ngrid,nslope) ! CO2 in the subslope [kg/m^2] |
|---|
| 164 | REAL, INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope |
|---|
| 165 | REAL, INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! boolean to check if there is flow in the mesh |
|---|
| 166 | ! Local |
|---|
| 167 | INTEGER ig,islope ! loop |
|---|
| 168 | REAL rho_co2(ngrid) ! density of CO2, temperature dependant [kg/m^3] |
|---|
| 169 | INTEGER iaval ! ice will be transfered here |
|---|
| 170 | |
|---|
| 171 | |
|---|
| 172 | |
|---|
| 173 | ! 0. Compute rho_co2 |
|---|
| 174 | DO ig = 1,ngrid |
|---|
| 175 | rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3 ! Mangan et al. 2017 |
|---|
| 176 | ENDDO |
|---|
| 177 | |
|---|
| 178 | ! 1. Compute the transfer of ice |
|---|
| 179 | |
|---|
| 180 | DO ig = 1,ngrid |
|---|
| 181 | DO islope = 1,nslope |
|---|
| 182 | IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground |
|---|
| 183 | ! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely on flat ground |
|---|
| 184 | IF(co2ice_slope(ig,islope).ge.rho_co2(ig)*hmax(ig,islope) * & |
|---|
| 185 | cos(pi*def_slope_mean(islope)/180.)) THEN |
|---|
| 186 | ! Second: determine the flatest slopes possible: |
|---|
| 187 | IF(islope.gt.iflat) THEN |
|---|
| 188 | iaval=islope-1 |
|---|
| 189 | ELSE |
|---|
| 190 | iaval=islope+1 |
|---|
| 191 | ENDIF |
|---|
| 192 | do while ((iaval.ne.iflat).and. & |
|---|
| 193 | (subslope_dist(ig,iaval).eq.0)) |
|---|
| 194 | IF(iaval.gt.iflat) THEN |
|---|
| 195 | iaval=iaval-1 |
|---|
| 196 | ELSE |
|---|
| 197 | iaval=iaval+1 |
|---|
| 198 | ENDIF |
|---|
| 199 | enddo |
|---|
| 200 | co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + & |
|---|
| 201 | (co2ice_slope(ig,islope) - rho_co2(ig)*hmax(ig,islope) * & |
|---|
| 202 | cos(pi*def_slope_mean(islope)/180.)) * & |
|---|
| 203 | subslope_dist(ig,islope)/subslope_dist(ig,iaval) * & |
|---|
| 204 | cos(pi*def_slope_mean(iaval)/180.) / & |
|---|
| 205 | cos(pi*def_slope_mean(islope)/180.) |
|---|
| 206 | |
|---|
| 207 | co2ice_slope(ig,islope)=rho_co2(ig)*hmax(ig,islope) * & |
|---|
| 208 | cos(pi*def_slope_mean(islope)/180.) |
|---|
| 209 | |
|---|
| 210 | flag_co2flow(ig,islope) = 1. |
|---|
| 211 | flag_co2flow_mesh(ig) = 1. |
|---|
| 212 | ENDIF ! co2ice > hmax |
|---|
| 213 | ENDIF ! iflat |
|---|
| 214 | ENDDO !islope |
|---|
| 215 | ENDDO !ig |
|---|
| 216 | RETURN |
|---|
| 217 | end subroutine |
|---|
| 218 | |
|---|
| 219 | subroutine computeTcond(timelen,ngrid,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) |
|---|
| 220 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 221 | !!! |
|---|
| 222 | !!! Purpose: Compute CO2 condensation temperature |
|---|
| 223 | !!! |
|---|
| 224 | !!! Author: LL |
|---|
| 225 | !!! |
|---|
| 226 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 227 | |
|---|
| 228 | use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2 |
|---|
| 229 | |
|---|
| 230 | implicit none |
|---|
| 231 | |
|---|
| 232 | ! arguments: |
|---|
| 233 | ! ---------- |
|---|
| 234 | |
|---|
| 235 | ! INPUT |
|---|
| 236 | INTEGER,INTENT(IN) :: timelen, ngrid ! # of timesample,physical points, subslopes |
|---|
| 237 | REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol] |
|---|
| 238 | REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa] |
|---|
| 239 | REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa] |
|---|
| 240 | REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration |
|---|
| 241 | ! OUTPUT |
|---|
| 242 | REAL,INTENT(OUT) :: Tcond(ngrid) ! Physical points : condensation temperature of CO2, yearly averaged |
|---|
| 243 | |
|---|
| 244 | ! LOCAL |
|---|
| 245 | |
|---|
| 246 | INTEGER :: ig,it ! for loop |
|---|
| 247 | REAL :: ave ! intermediate to compute average |
|---|
| 248 | |
|---|
| 249 | !!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 250 | |
|---|
| 251 | |
|---|
| 252 | DO ig = 1,ngrid |
|---|
| 253 | ave = 0 |
|---|
| 254 | DO it = 1,timelen |
|---|
| 255 | 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)) |
|---|
| 256 | ENDDO |
|---|
| 257 | Tcond(ig) = ave/timelen |
|---|
| 258 | ENDDO |
|---|
| 259 | RETURN |
|---|
| 260 | |
|---|
| 261 | |
|---|
| 262 | end subroutine |
|---|
| 263 | end module |
|---|
