| 1 | subroutine spreadglacier_paleo(ngrid,nq,pqsurf, & |
|---|
| 2 | phisfi0,timstep,ptsurf) |
|---|
| 3 | |
|---|
| 4 | use callkeys_mod, only: carbox, methane |
|---|
| 5 | use comcstfi_mod, only: g |
|---|
| 6 | use comgeomfi_h |
|---|
| 7 | use geometry_mod, only: latitude, longitude, cell_area |
|---|
| 8 | use tracer_h, only: igcm_n2,igcm_ch4_ice,igcm_co_ice |
|---|
| 9 | |
|---|
| 10 | implicit none |
|---|
| 11 | |
|---|
| 12 | !================================================================== |
|---|
| 13 | ! Purpose |
|---|
| 14 | ! ------- |
|---|
| 15 | ! Spreading the glacier : N2, cf Umurhan et al 2017 |
|---|
| 16 | ! |
|---|
| 17 | ! Inputs |
|---|
| 18 | ! ------ |
|---|
| 19 | ! ngrid Number of vertical columns |
|---|
| 20 | ! pqsurf(ngrid,nq) N2 ice tracer on surface kg/m2 |
|---|
| 21 | ! phisfi0 = phisfinew the geopotential of the bedrock |
|---|
| 22 | ! below the N2 ice |
|---|
| 23 | ! ptsurf surface temperature |
|---|
| 24 | ! timstep dstep = timestep in pluto day for the call of glacial flow |
|---|
| 25 | |
|---|
| 26 | ! Outputs |
|---|
| 27 | ! ------- |
|---|
| 28 | ! pqsurf(ngrid,nq) new value for N2 ice tracer on surface kg/m2 |
|---|
| 29 | ! |
|---|
| 30 | ! Authors |
|---|
| 31 | ! ------- |
|---|
| 32 | ! Tanguy Bertrand (2015,2022) |
|---|
| 33 | ! |
|---|
| 34 | !================================================================== |
|---|
| 35 | |
|---|
| 36 | #include "dimensions.h" |
|---|
| 37 | |
|---|
| 38 | !----------------------------------------------------------------------- |
|---|
| 39 | ! Arguments |
|---|
| 40 | |
|---|
| 41 | INTEGER ngrid, nq, ig, i |
|---|
| 42 | REAL :: pqsurf(ngrid,nq) |
|---|
| 43 | REAL :: phisfi0(ngrid) |
|---|
| 44 | REAL :: ptsurf(ngrid) |
|---|
| 45 | REAL timstep |
|---|
| 46 | !----------------------------------------------------------------------- |
|---|
| 47 | ! Local variables |
|---|
| 48 | REAL tgla(ngrid),tbase(ngrid) !temperature at the base of glacier different of ptsurf |
|---|
| 49 | REAL :: zdqsurf(ngrid) ! tendency of qsurf |
|---|
| 50 | REAL massflow ! function |
|---|
| 51 | REAL hlim,qlim,difflim,diff,stock,H0,totmasstmp |
|---|
| 52 | INTEGER compt !compteur |
|---|
| 53 | INTEGER slid ! option slid 0 or 1 |
|---|
| 54 | INTEGER :: edge(4) ! index of the adjecent points, N,S,E,W |
|---|
| 55 | INTEGER cas,inddeb,indfin |
|---|
| 56 | REAL secu,dqch4,dqco |
|---|
| 57 | REAL :: masstmp(ngrid) |
|---|
| 58 | |
|---|
| 59 | REAL :: masstot |
|---|
| 60 | |
|---|
| 61 | !---------------- INPUT ------------------------------------------------ |
|---|
| 62 | |
|---|
| 63 | difflim=0.5 ! limit height (m) between two adjacent point to start spreading the glacier |
|---|
| 64 | zdqsurf(:)=0. |
|---|
| 65 | !--------------- Dimensions ------------------------------------- |
|---|
| 66 | |
|---|
| 67 | ! distance between two adjacent latitude |
|---|
| 68 | !distlim_lat=pi*rad/jjm/1000. ! km |
|---|
| 69 | !distlim_diag=(distlim_lat*distlim_lat+distlim_lon*distlim_lon)**0.5 |
|---|
| 70 | |
|---|
| 71 | !! Threshold for ice thickness for which we activate the glacial flow |
|---|
| 72 | hlim=10. !m |
|---|
| 73 | qlim=hlim*1000. ! kg |
|---|
| 74 | !! Security for not depleted all ice too fast in one timestep |
|---|
| 75 | secu=4 |
|---|
| 76 | |
|---|
| 77 | !************************************* |
|---|
| 78 | ! Loop over all points |
|---|
| 79 | !************************************* |
|---|
| 80 | DO ig=1,ngrid |
|---|
| 81 | !if (ig.eq.ngrid) then |
|---|
| 82 | ! print*, 'qpole=',pqsurf(ig,igcm_n2),qlim |
|---|
| 83 | !endif |
|---|
| 84 | |
|---|
| 85 | !************************************* |
|---|
| 86 | ! if significant amount of ice, where qsurf > qlim |
|---|
| 87 | !************************************* |
|---|
| 88 | IF (pqsurf(ig,igcm_n2).gt.qlim) THEN |
|---|
| 89 | |
|---|
| 90 | !************************************* |
|---|
| 91 | ! temp glacier with gradient 15 K / km |
|---|
| 92 | !************************************* |
|---|
| 93 | ! surface temperature pts |
|---|
| 94 | tgla(ig)=ptsurf(ig) |
|---|
| 95 | ! temperature at the base (we neglect convection) |
|---|
| 96 | tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al. |
|---|
| 97 | if (tbase(ig).gt.63.) then |
|---|
| 98 | slid=1 ! activate slide |
|---|
| 99 | else |
|---|
| 100 | slid=0 |
|---|
| 101 | endif |
|---|
| 102 | |
|---|
| 103 | !************************************* |
|---|
| 104 | ! Selection of the adjacent index depending on the grid point |
|---|
| 105 | !************************************* |
|---|
| 106 | !! poles |
|---|
| 107 | IF (ig.eq.1) then |
|---|
| 108 | cas=0 |
|---|
| 109 | inddeb=2 ! First adj grid point |
|---|
| 110 | indfin=iim+1 ! Last adj grid point |
|---|
| 111 | ELSEIF (ig.eq.ngrid) then |
|---|
| 112 | cas=10 |
|---|
| 113 | inddeb=ngrid-iim |
|---|
| 114 | indfin=ngrid-1 |
|---|
| 115 | !! 9 other cases: edges East, west, or in the center , at edges north south or in the center |
|---|
| 116 | ELSEIF (mod(ig,iim).eq.2) then ! Edge east : N,S,W,E |
|---|
| 117 | IF (ig.eq.2) then |
|---|
| 118 | cas=1 |
|---|
| 119 | edge = (/1, ig+iim,ig-1+iim,ig+1 /) |
|---|
| 120 | ELSEIF (ig.eq.ngrid-iim) then |
|---|
| 121 | cas=3 |
|---|
| 122 | edge = (/ig-iim, ngrid,ig-1+iim,ig+1 /) |
|---|
| 123 | ELSE |
|---|
| 124 | cas=2 |
|---|
| 125 | edge = (/ig-iim, ig+iim,ig-1+iim,ig+1 /) |
|---|
| 126 | ENDIF |
|---|
| 127 | ELSEIF (mod(ig,iim).eq.1) then ! Edge west |
|---|
| 128 | IF (ig.eq.iim+1) then |
|---|
| 129 | cas=7 |
|---|
| 130 | edge = (/1,ig+iim,ig-1,ig+1-iim /) |
|---|
| 131 | ELSEIF (ig.eq.ngrid-1) then |
|---|
| 132 | cas=9 |
|---|
| 133 | edge = (/ig-iim,ngrid,ig-1,ig+1-iim /) |
|---|
| 134 | ELSE |
|---|
| 135 | cas=8 |
|---|
| 136 | edge = (/ig-iim,ig+iim,ig-1,ig+1-iim /) |
|---|
| 137 | ENDIF |
|---|
| 138 | ELSE |
|---|
| 139 | IF ((ig.gt.2).and.(ig.lt.iim+1)) then |
|---|
| 140 | cas=4 |
|---|
| 141 | edge = (/1,ig+iim,ig-1,ig+1 /) |
|---|
| 142 | ELSEIF ((ig.gt.ngrid-iim).and.(ig.lt.ngrid-1)) then |
|---|
| 143 | cas=6 |
|---|
| 144 | edge = (/ig-iim,ngrid,ig-1,ig+1 /) |
|---|
| 145 | ELSE |
|---|
| 146 | cas=5 |
|---|
| 147 | edge = (/ig-iim,ig+iim,ig-1,ig+1 /) |
|---|
| 148 | ENDIF |
|---|
| 149 | ENDIF |
|---|
| 150 | |
|---|
| 151 | masstmp(:)=0. ! mass to be transferred |
|---|
| 152 | totmasstmp=0. ! total mass to be transferred |
|---|
| 153 | H0=phisfi0(ig)/g+pqsurf(ig,igcm_n2)/1000. ! height of ice at ig (m) |
|---|
| 154 | |
|---|
| 155 | !************************************* |
|---|
| 156 | !!!!!!!!!!!!!!!!! POLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 157 | !************************************* |
|---|
| 158 | IF ((cas.eq.0).or.(cas.eq.10)) then ! Mean over all latitudes near the pole |
|---|
| 159 | |
|---|
| 160 | !call testconservmass2d(ngrid,pqsurf(:,1)) |
|---|
| 161 | !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point |
|---|
| 162 | DO i=inddeb,indfin |
|---|
| 163 | diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! Height difference between ig and adjacent points (m) |
|---|
| 164 | if (diff.gt.difflim) then |
|---|
| 165 | masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s) |
|---|
| 166 | else |
|---|
| 167 | masstmp(i)=0. |
|---|
| 168 | endif |
|---|
| 169 | totmasstmp=totmasstmp+masstmp(i) ! mass totale to be transferred if assume only one adjecent point |
|---|
| 170 | ENDDO |
|---|
| 171 | if (totmasstmp.gt.0.) then |
|---|
| 172 | !! 2) Available mass to be transferred |
|---|
| 173 | stock=maxval(masstmp(:))/secu ! kg/m2/s |
|---|
| 174 | !! 3) Real amounts of ice to be transferred : |
|---|
| 175 | masstmp(:)=masstmp(:)/totmasstmp*stock*cell_area(ig)/cell_area(inddeb) !kg/m2/s all area around the pole are the same |
|---|
| 176 | !! 4) Tendencies |
|---|
| 177 | zdqsurf(ig)=-stock !kg/m2/s |
|---|
| 178 | !! 5) Update PQSURF |
|---|
| 179 | |
|---|
| 180 | ! move CH4 and CO too |
|---|
| 181 | if (methane) then |
|---|
| 182 | !! Variation of CH4 ice |
|---|
| 183 | dqch4=pqsurf(ig,igcm_ch4_ice)*(zdqsurf(ig)/pqsurf(ig,igcm_n2)*timstep) ! amout of ch4 (kg/m2) to be removed at grid ig (negative) |
|---|
| 184 | !! remove CH4 |
|---|
| 185 | pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 |
|---|
| 186 | !! add CH4 |
|---|
| 187 | do i=inddeb,indfin |
|---|
| 188 | pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 |
|---|
| 189 | enddo |
|---|
| 190 | endif |
|---|
| 191 | if (carbox) then |
|---|
| 192 | !! Variation of CO ice |
|---|
| 193 | dqco=pqsurf(ig,igcm_co_ice)*(zdqsurf(ig)/pqsurf(ig,igcm_n2)*timstep) ! amout of co (kg/m2) to be removed at grid ig (negative) |
|---|
| 194 | !! remove CO |
|---|
| 195 | pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco |
|---|
| 196 | !! add CO |
|---|
| 197 | do i=inddeb,indfin |
|---|
| 198 | pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco |
|---|
| 199 | enddo |
|---|
| 200 | endif |
|---|
| 201 | |
|---|
| 202 | ! Add N2 |
|---|
| 203 | do i=inddeb,indfin |
|---|
| 204 | pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep |
|---|
| 205 | enddo |
|---|
| 206 | ! Remove N2 |
|---|
| 207 | pqsurf(ig,igcm_n2)=pqsurf(ig,igcm_n2)+zdqsurf(ig)*timstep |
|---|
| 208 | |
|---|
| 209 | endif |
|---|
| 210 | !call testconservmass2d(ngrid,pqsurf(:,1)) |
|---|
| 211 | |
|---|
| 212 | !!!!!!!!!!!!!!!!!!!! NOT THE POLE !!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 213 | ! here: ig = grid point with large amount of ice |
|---|
| 214 | ! Loop on each adjacent point |
|---|
| 215 | ! looking for adjacent point |
|---|
| 216 | ELSE |
|---|
| 217 | |
|---|
| 218 | !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point |
|---|
| 219 | DO compt=1,4 |
|---|
| 220 | i=edge(compt) |
|---|
| 221 | diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! (m) |
|---|
| 222 | if (diff.gt.difflim) then |
|---|
| 223 | masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) |
|---|
| 224 | else |
|---|
| 225 | masstmp(i)=0. |
|---|
| 226 | endif |
|---|
| 227 | totmasstmp=totmasstmp+masstmp(i) |
|---|
| 228 | ENDDO |
|---|
| 229 | if (totmasstmp.gt.0) then |
|---|
| 230 | !! 2) Available mass to be transferred |
|---|
| 231 | stock=maxval(masstmp(:))/secu ! kg/m2/s |
|---|
| 232 | !! 3) Real amounts of ice to be transferred : |
|---|
| 233 | DO compt=1,4 |
|---|
| 234 | i=edge(compt) |
|---|
| 235 | masstmp(i)=masstmp(i)/totmasstmp*stock*cell_area(ig)/cell_area(i) !kg/m2/s |
|---|
| 236 | ENDDO |
|---|
| 237 | !! 4) Tendencies |
|---|
| 238 | zdqsurf(ig)=-stock |
|---|
| 239 | !! 5) Update PQSURF |
|---|
| 240 | |
|---|
| 241 | ! move CH4 and CO too |
|---|
| 242 | if (methane) then |
|---|
| 243 | !! Variation of CH4 ice |
|---|
| 244 | dqch4=pqsurf(ig,igcm_ch4_ice)*(zdqsurf(ig)/pqsurf(ig,igcm_n2)*timstep) ! amout of ch4 (kg/m2) to be removed at grid ig (negative) |
|---|
| 245 | !! remove CH4 |
|---|
| 246 | pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 |
|---|
| 247 | !! add CH4 |
|---|
| 248 | DO compt=1,4 |
|---|
| 249 | i=edge(compt) |
|---|
| 250 | pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 |
|---|
| 251 | enddo |
|---|
| 252 | endif |
|---|
| 253 | if (carbox) then |
|---|
| 254 | !! Variation of CO ice |
|---|
| 255 | dqco=pqsurf(ig,igcm_co_ice)*(zdqsurf(ig)/pqsurf(ig,igcm_n2)*timstep) ! amout of co (kg/m2) to be removed at grid ig (negative) |
|---|
| 256 | !! remove CO |
|---|
| 257 | pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco |
|---|
| 258 | !! add CO |
|---|
| 259 | DO compt=1,4 |
|---|
| 260 | i=edge(compt) |
|---|
| 261 | pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco |
|---|
| 262 | enddo |
|---|
| 263 | endif |
|---|
| 264 | |
|---|
| 265 | ! Add N2 |
|---|
| 266 | totmasstmp=0. |
|---|
| 267 | DO compt=1,4 |
|---|
| 268 | i=edge(compt) |
|---|
| 269 | pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep |
|---|
| 270 | totmasstmp=totmasstmp+masstmp(i)*cell_area(i) !! TB22 tmp |
|---|
| 271 | enddo |
|---|
| 272 | ! Remove N2 |
|---|
| 273 | pqsurf(ig,igcm_n2)=pqsurf(ig,igcm_n2)+zdqsurf(ig)*timstep |
|---|
| 274 | |
|---|
| 275 | endif |
|---|
| 276 | |
|---|
| 277 | if (pqsurf(ig,igcm_n2).lt.0) then |
|---|
| 278 | print*, pqsurf(ig,igcm_n2) |
|---|
| 279 | write(*,*) 'bug in spreadglacier' |
|---|
| 280 | stop |
|---|
| 281 | endif |
|---|
| 282 | |
|---|
| 283 | ENDIF ! cas |
|---|
| 284 | ENDIF ! qlim |
|---|
| 285 | |
|---|
| 286 | ENDDO ! ig |
|---|
| 287 | |
|---|
| 288 | end subroutine spreadglacier_paleo |
|---|
| 289 | |
|---|
| 290 | real function dist_pluto(lat1,lon1,lat2,lon2) |
|---|
| 291 | use comcstfi_mod, only: rad |
|---|
| 292 | implicit none |
|---|
| 293 | real, intent(in) :: lat1,lon1,lat2,lon2 |
|---|
| 294 | real dlon,dlat |
|---|
| 295 | real a,c |
|---|
| 296 | real radi |
|---|
| 297 | radi=rad/1000. |
|---|
| 298 | |
|---|
| 299 | dlon = lon2 - lon1 |
|---|
| 300 | dlat = lat2 - lat1 |
|---|
| 301 | a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 |
|---|
| 302 | c = 2 * atan2( sqrt(a), sqrt(1-a) ) |
|---|
| 303 | dist_pluto = radi * c !! km |
|---|
| 304 | return |
|---|
| 305 | end function dist_pluto |
|---|
| 306 | |
|---|
| 307 | real function massflow(ig1,ig2,pts,dif,pqs,dt,slid) |
|---|
| 308 | !! Mass of ice to be transferred from one grid point to another |
|---|
| 309 | use comgeomfi_h |
|---|
| 310 | use geometry_mod, only: latitude, longitude, cell_area |
|---|
| 311 | use comcstfi_mod, only: g, rad |
|---|
| 312 | use tracer_h, only: rho_n2 |
|---|
| 313 | |
|---|
| 314 | implicit none |
|---|
| 315 | #include "dimensions.h" |
|---|
| 316 | |
|---|
| 317 | real, intent(in) :: pts,dif,pqs,dt |
|---|
| 318 | integer, intent(in) :: ig1,ig2,slid |
|---|
| 319 | real lref,ac,n,theta,hdelta,ha,tau,qdry,qsl,usl0 |
|---|
| 320 | REAL dist_pluto ! function |
|---|
| 321 | |
|---|
| 322 | lref=dist_pluto(latitude(ig2),longitude(ig2),latitude(ig1),longitude(ig1))*1000. ! m |
|---|
| 323 | usl0=6.3e-6 ! m/s |
|---|
| 324 | ac=0.005*exp(9.37-422./pts) !0.005*exp(422./45.-422./pts) |
|---|
| 325 | n=2.1+0.0155*(pts-45.) |
|---|
| 326 | theta=atan(dif/(lref)) |
|---|
| 327 | hdelta=pts/20. |
|---|
| 328 | ha=pts*pts/8440. !pts*pts/20./422. |
|---|
| 329 | qdry=ac*(rho_n2*g*1.e-6)**n*(pqs/1000.)**(n+2)/(n+2)*tan(theta)**(n-1)/(1+tan(theta)*tan(theta))**(n/2.) |
|---|
| 330 | qdry=qdry*0.5*exp( (pqs/1000./ha) / ( 1+ pqs /hdelta) ) |
|---|
| 331 | qsl=(pqs/1000.)/abs(tan(theta))*usl0*slid |
|---|
| 332 | tau=pqs/1000.*lref/abs((qdry+qsl)*tan(theta)) |
|---|
| 333 | massflow=pqs*(1.-exp(-dt/tau))/dt |
|---|
| 334 | return |
|---|
| 335 | end function massflow |
|---|
| 336 | |
|---|
| 337 | |
|---|
| 338 | |
|---|
| 339 | |
|---|
