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