[2995] | 1 | module glaciers_mod |
---|
| 2 | |
---|
| 3 | implicit none |
---|
| 4 | LOGICAL co2glaciersflow ! True by default, to compute co2 ice flow. Read in pem.def |
---|
| 5 | LOGICAL h2oglaciersflow ! True by default, to compute co2 ice flow. Read in pem.def |
---|
| 6 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 7 | !!! |
---|
| 8 | !!! Purpose: Compute CO2 glacier flows |
---|
| 9 | !!! |
---|
| 10 | !!! Author: LL |
---|
| 11 | !!! |
---|
| 12 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 13 | |
---|
| 14 | contains |
---|
| 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,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 | IMPLICIT NONE |
---|
| 29 | |
---|
| 30 | ! arguments |
---|
| 31 | ! --------- |
---|
| 32 | |
---|
| 33 | ! Inputs: |
---|
| 34 | INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! number of time sample, physical points, subslopes, index of the flat subslope |
---|
| 35 | 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 |
---|
| 36 | REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol] |
---|
| 37 | REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical x Time field: surface pressure given by the GCM [Pa] |
---|
| 38 | REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surface pressure from the GCM [Pa] |
---|
| 39 | REAL,INTENT(IN) :: global_ave_ps_PEM ! global averaged surface pressure during the PEM iteration [Pa] |
---|
| 40 | |
---|
| 41 | ! Ouputs: |
---|
| 42 | REAL,INTENT(INOUT) :: co2ice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] |
---|
| 43 | REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes |
---|
| 44 | REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid) ! same but within the mesh |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | ! Local |
---|
| 48 | REAL :: Tcond(ngrid,nslope) ! Physical field: CO2 condensation temperature [K] |
---|
| 49 | REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2 glacier before flow |
---|
| 50 | |
---|
| 51 | !----------------------------- |
---|
| 52 | call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) |
---|
| 53 | |
---|
| 54 | call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax) |
---|
| 55 | |
---|
| 56 | call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh) |
---|
[3002] | 57 | RETURN |
---|
[2995] | 58 | end subroutine |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | |
---|
| 63 | |
---|
| 64 | subroutine h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh) |
---|
| 65 | |
---|
| 66 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 67 | !!! |
---|
| 68 | !!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do |
---|
| 69 | !!! the ice transfer |
---|
| 70 | !!! |
---|
| 71 | !!! |
---|
| 72 | !!! Author: LL |
---|
| 73 | !!! |
---|
| 74 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | IMPLICIT NONE |
---|
| 78 | |
---|
| 79 | ! arguments |
---|
| 80 | ! --------- |
---|
| 81 | |
---|
| 82 | ! Inputs: |
---|
| 83 | INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! number of time sample, physical points, subslopes, index of the flat subslope |
---|
| 84 | 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 |
---|
| 85 | REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Ice Temperature [K] |
---|
| 86 | ! Ouputs: |
---|
| 87 | REAL,INTENT(INOUT) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2] |
---|
| 88 | REAL,INTENT(INOUT) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes |
---|
| 89 | REAL,INTENT(INOUT) :: flag_h2oflow_mesh(ngrid) ! same but within the mesh |
---|
| 90 | ! Local |
---|
| 91 | REAL :: hmax(ngrid,nslope) ! Physical x Slope field: maximum thickness for co2 glacier before flow |
---|
| 92 | |
---|
| 93 | !----------------------------- |
---|
| 94 | |
---|
| 95 | call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax) |
---|
| 96 | call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh) |
---|
| 97 | |
---|
| 98 | RETURN |
---|
| 99 | end subroutine |
---|
| 100 | |
---|
| 101 | |
---|
| 102 | subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax) |
---|
| 103 | |
---|
[3076] | 104 | use abort_pem_mod, only: abort_pem |
---|
[3082] | 105 | #ifndef CPP_STD |
---|
| 106 | use comcstfi_h, only: pi, g |
---|
| 107 | #else |
---|
| 108 | use comcstfi_mod, only: pi, g |
---|
| 109 | #endif |
---|
[2995] | 110 | |
---|
| 111 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 112 | !!! |
---|
| 113 | !!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle |
---|
| 114 | !!! before initating flow |
---|
| 115 | !!! |
---|
| 116 | !!! Author: LL,based on work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022) |
---|
| 117 | !!! |
---|
| 118 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 119 | |
---|
| 120 | IMPLICIT NONE |
---|
| 121 | |
---|
| 122 | ! arguments |
---|
| 123 | ! -------- |
---|
| 124 | |
---|
| 125 | ! Inputs |
---|
| 126 | INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes |
---|
| 127 | INTEGER,INTENT(IN) :: iflat ! index of the flat subslope |
---|
| 128 | REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg] |
---|
| 129 | REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Physical field: ice temperature [K] |
---|
[3002] | 130 | character(len=3), INTENT(IN) :: name_ice ! Nature of the ice |
---|
[2995] | 131 | ! Outputs |
---|
| 132 | REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum thickness before flaw [m] |
---|
| 133 | ! Local |
---|
| 134 | 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 |
---|
| 135 | REAL :: rho(ngrid,nslope) ! co2 ice density [kg/m^3] |
---|
| 136 | INTEGER :: ig,islope ! loop variables |
---|
| 137 | REAL :: slo_angle |
---|
| 138 | |
---|
| 139 | ! 1. Compute rho |
---|
| 140 | if(name_ice.eq."co2") then |
---|
| 141 | DO ig = 1,ngrid |
---|
| 142 | DO islope = 1,nslope |
---|
| 143 | rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3 ! Mangan et al. 2017 |
---|
| 144 | tau_d = 5.e3 |
---|
| 145 | ENDDO |
---|
| 146 | ENDDO |
---|
| 147 | elseif (name_ice.eq."h2o") then |
---|
| 148 | DO ig = 1,ngrid |
---|
| 149 | DO islope = 1,nslope |
---|
| 150 | rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 |
---|
| 151 | tau_d = 1.e5 |
---|
| 152 | ENDDO |
---|
| 153 | ENDDO |
---|
| 154 | else |
---|
| 155 | call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1) |
---|
| 156 | endif |
---|
| 157 | |
---|
| 158 | ! 3. Compute max thickness |
---|
| 159 | DO ig = 1,ngrid |
---|
| 160 | DO islope = 1,nslope |
---|
| 161 | if(islope.eq.iflat) then |
---|
| 162 | hmax(ig,islope) = 1.e8 |
---|
| 163 | else |
---|
| 164 | slo_angle = abs(def_slope_mean(islope)*pi/180.) |
---|
| 165 | hmax(ig,islope) = tau_d/(rho(ig,islope)*g*slo_angle) |
---|
| 166 | endif |
---|
| 167 | ENDDO |
---|
| 168 | ENDDO |
---|
| 169 | RETURN |
---|
| 170 | |
---|
| 171 | end subroutine |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | |
---|
| 175 | |
---|
| 176 | subroutine transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh) |
---|
| 177 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 178 | !!! |
---|
| 179 | !!! Purpose: Transfer the excess of ice from one subslope to another |
---|
| 180 | !!! No transfer between mesh at the time |
---|
| 181 | !!! Author: LL |
---|
| 182 | !!! |
---|
| 183 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 184 | |
---|
[3076] | 185 | use abort_pem_mod, only: abort_pem |
---|
[3082] | 186 | #ifndef CPP_STD |
---|
| 187 | use comcstfi_h, only: pi |
---|
| 188 | #else |
---|
| 189 | use comcstfi_mod, only: pi |
---|
| 190 | #endif |
---|
[2995] | 191 | |
---|
| 192 | implicit none |
---|
| 193 | |
---|
| 194 | ! arguments |
---|
| 195 | ! -------- |
---|
| 196 | |
---|
| 197 | ! Inputs |
---|
| 198 | INTEGER, INTENT(IN) :: ngrid,nslope !# of physical points and subslope |
---|
| 199 | INTEGER, INTENT(IN) :: iflat ! index of the flat subslope |
---|
| 200 | REAL, INTENT(IN) :: subslope_dist(ngrid,nslope) ! Distribution of the subgrid slopes within the mesh |
---|
| 201 | REAL, INTENT(IN) :: def_slope_mean(nslope) ! values of the subgrid slopes |
---|
| 202 | REAL, INTENT(IN) :: hmax(ngrid,nslope) ! maximum height of the glaciers before initiating flow [m] |
---|
| 203 | REAL, INTENT(IN) :: Tice(ngrid,nslope) ! Ice temperature[K] |
---|
[3002] | 204 | character(len=3), INTENT(IN) :: name_ice ! Nature of the ice |
---|
[2995] | 205 | |
---|
| 206 | ! Outputs |
---|
| 207 | REAL, INTENT(INOUT) :: qice(ngrid,nslope) ! CO2 in the subslope [kg/m^2] |
---|
| 208 | REAL, INTENT(INOUT) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope |
---|
| 209 | REAL, INTENT(INOUT) :: flag_flowmesh(ngrid) ! boolean to check if there is flow in the mesh |
---|
| 210 | ! Local |
---|
| 211 | INTEGER ig,islope ! loop |
---|
| 212 | REAL rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3] |
---|
| 213 | INTEGER iaval ! ice will be transfered here |
---|
| 214 | |
---|
| 215 | ! 0. Compute rho |
---|
| 216 | if(name_ice.eq."co2") then |
---|
| 217 | DO ig = 1,ngrid |
---|
| 218 | DO islope = 1,nslope |
---|
| 219 | rho(ig,islope) = (1.72391 - 2.53e-4*Tice(ig,islope)-2.87*1e-7*Tice(ig,islope)**2)*1e3 ! Mangan et al. 2017 |
---|
| 220 | ENDDO |
---|
| 221 | ENDDO |
---|
| 222 | elseif (name_ice.eq."h2o") then |
---|
| 223 | DO ig = 1,ngrid |
---|
| 224 | DO islope = 1,nslope |
---|
| 225 | rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 |
---|
| 226 | ENDDO |
---|
| 227 | ENDDO |
---|
| 228 | else |
---|
| 229 | call abort_pem("PEM - Transfer ice","Name of ice is not co2 or h2o",1) |
---|
| 230 | endif |
---|
| 231 | |
---|
| 232 | ! 1. Compute the transfer of ice |
---|
| 233 | |
---|
| 234 | DO ig = 1,ngrid |
---|
| 235 | DO islope = 1,nslope |
---|
| 236 | IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground |
---|
| 237 | ! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate infinitely on flat ground |
---|
| 238 | IF(qice(ig,islope).ge.rho(ig,islope)*hmax(ig,islope) * & |
---|
| 239 | cos(pi*def_slope_mean(islope)/180.)) THEN |
---|
| 240 | ! Second: determine the flatest slopes possible: |
---|
| 241 | IF(islope.gt.iflat) THEN |
---|
| 242 | iaval=islope-1 |
---|
| 243 | ELSE |
---|
| 244 | iaval=islope+1 |
---|
| 245 | ENDIF |
---|
| 246 | do while ((iaval.ne.iflat).and. & |
---|
| 247 | (subslope_dist(ig,iaval).eq.0)) |
---|
| 248 | IF(iaval.gt.iflat) THEN |
---|
| 249 | iaval=iaval-1 |
---|
| 250 | ELSE |
---|
| 251 | iaval=iaval+1 |
---|
| 252 | ENDIF |
---|
| 253 | enddo |
---|
| 254 | qice(ig,iaval) = qice(ig,iaval) + & |
---|
| 255 | (qice(ig,islope) - rho(ig,islope)*hmax(ig,islope) * & |
---|
| 256 | cos(pi*def_slope_mean(islope)/180.)) * & |
---|
| 257 | subslope_dist(ig,islope)/subslope_dist(ig,iaval) * & |
---|
| 258 | cos(pi*def_slope_mean(iaval)/180.) / & |
---|
| 259 | cos(pi*def_slope_mean(islope)/180.) |
---|
| 260 | |
---|
| 261 | qice(ig,islope)=rho(ig,islope)*hmax(ig,islope) * & |
---|
| 262 | cos(pi*def_slope_mean(islope)/180.) |
---|
| 263 | |
---|
| 264 | flag_flow(ig,islope) = 1. |
---|
| 265 | flag_flowmesh(ig) = 1. |
---|
| 266 | ENDIF ! co2ice > hmax |
---|
| 267 | ENDIF ! iflat |
---|
| 268 | ENDDO !islope |
---|
| 269 | ENDDO !ig |
---|
| 270 | RETURN |
---|
| 271 | end subroutine |
---|
| 272 | |
---|
| 273 | subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond) |
---|
| 274 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 275 | !!! |
---|
| 276 | !!! Purpose: Compute CO2 condensation temperature |
---|
| 277 | !!! |
---|
| 278 | !!! Author: LL |
---|
| 279 | !!! |
---|
| 280 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 281 | |
---|
| 282 | use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2 |
---|
| 283 | |
---|
| 284 | implicit none |
---|
| 285 | |
---|
| 286 | ! arguments: |
---|
| 287 | ! ---------- |
---|
| 288 | |
---|
| 289 | ! INPUT |
---|
| 290 | INTEGER,INTENT(IN) :: timelen, ngrid,nslope ! # of timesample,physical points, subslopes |
---|
| 291 | REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol] |
---|
| 292 | REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa] |
---|
| 293 | REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa] |
---|
| 294 | REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration |
---|
| 295 | ! OUTPUT |
---|
| 296 | REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged |
---|
| 297 | |
---|
| 298 | ! LOCAL |
---|
| 299 | |
---|
| 300 | INTEGER :: ig,it,islope ! for loop |
---|
| 301 | REAL :: ave ! intermediate to compute average |
---|
| 302 | |
---|
| 303 | !!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 304 | |
---|
| 305 | |
---|
| 306 | DO ig = 1,ngrid |
---|
| 307 | ave = 0 |
---|
| 308 | DO it = 1,timelen |
---|
| 309 | 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)) |
---|
| 310 | ENDDO |
---|
| 311 | DO islope = 1,nslope |
---|
| 312 | Tcond(ig,islope) = ave/timelen |
---|
| 313 | ENDDO |
---|
| 314 | ENDDO |
---|
| 315 | RETURN |
---|
| 316 | |
---|
| 317 | |
---|
| 318 | end subroutine |
---|
| 319 | end module |
---|