| 1 | subroutine spreadglacier_paleo(ngrid,nq,pqsurf, & |
|---|
| 2 | phisfi0,timstep,ptsurf) |
|---|
| 3 | |
|---|
| 4 | use radinc_h, only : naerkind |
|---|
| 5 | use comgeomfi_h |
|---|
| 6 | implicit none |
|---|
| 7 | |
|---|
| 8 | !================================================================== |
|---|
| 9 | ! Purpose |
|---|
| 10 | ! ------- |
|---|
| 11 | ! Spreading the glacier : N2, cf Umurhan et al 2017 |
|---|
| 12 | ! |
|---|
| 13 | ! Inputs |
|---|
| 14 | ! ------ |
|---|
| 15 | ! ngrid Number of vertical columns |
|---|
| 16 | ! pqsurf(ngrid,nq) N2 ice tracer on surface kg/m2 |
|---|
| 17 | ! phisfi0 = phisfinew the geopotential of the bedrock |
|---|
| 18 | ! below the N2 ice |
|---|
| 19 | ! ptsurf surface temperature |
|---|
| 20 | ! timstep dstep = timestep in pluto day for the call of glacial flow |
|---|
| 21 | |
|---|
| 22 | ! Outputs |
|---|
| 23 | ! ------- |
|---|
| 24 | ! pqsurf(ngrid,nq) new value for N2 ice tracer on surface kg/m2 |
|---|
| 25 | ! |
|---|
| 26 | ! Authors |
|---|
| 27 | ! ------- |
|---|
| 28 | ! Tanguy Bertrand (2015) |
|---|
| 29 | ! |
|---|
| 30 | !================================================================== |
|---|
| 31 | |
|---|
| 32 | #include "dimensions.h" |
|---|
| 33 | #include "dimphys.h" |
|---|
| 34 | #include "comcstfi.h" |
|---|
| 35 | #include "surfdat.h" |
|---|
| 36 | #include "comvert.h" |
|---|
| 37 | #include "callkeys.h" |
|---|
| 38 | #include "tracer.h" |
|---|
| 39 | |
|---|
| 40 | !----------------------------------------------------------------------- |
|---|
| 41 | ! Arguments |
|---|
| 42 | |
|---|
| 43 | INTEGER ngrid, nq, ig, i |
|---|
| 44 | REAL :: pqsurf(ngrid,nq) |
|---|
| 45 | REAL :: phisfi0(ngrid) |
|---|
| 46 | REAL :: ptsurf(ngrid) |
|---|
| 47 | REAL timstep |
|---|
| 48 | !----------------------------------------------------------------------- |
|---|
| 49 | ! Local variables |
|---|
| 50 | REAL tgla(ngrid),tbase(ngrid) !temperature at the base of glacier different of ptsurf |
|---|
| 51 | REAL :: zdqsurf(ngrid) ! tendency of qsurf |
|---|
| 52 | REAL massflow ! function |
|---|
| 53 | REAL hlim,qlim,difflim,diff,stock,H0,totmasstmp |
|---|
| 54 | INTEGER compt !compteur |
|---|
| 55 | INTEGER slid ! option slid 0 or 1 |
|---|
| 56 | INTEGER :: edge(4) ! index of the adjecent points, N,S,E,W |
|---|
| 57 | INTEGER cas,inddeb,indfin |
|---|
| 58 | REAL secu,dqch4,dqco |
|---|
| 59 | REAL :: masstmp(ngrid) |
|---|
| 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 |
|---|
| 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 | !print*, 'APPEL cas', cas |
|---|
| 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 | !print*, 'APPEL1',cas, H0 |
|---|
| 162 | !call testconservmass2d(ngrid,pqsurf(:,1)) |
|---|
| 163 | !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point |
|---|
| 164 | DO i=inddeb,indfin |
|---|
| 165 | !print*, 'area i=',area(i),area(ig),i,ig |
|---|
| 166 | diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! Height difference between ig and adjacent points (m) |
|---|
| 167 | !print*, 'diff=',diff |
|---|
| 168 | if (diff.gt.difflim) then |
|---|
| 169 | masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s) |
|---|
| 170 | else |
|---|
| 171 | masstmp(i)=0. |
|---|
| 172 | endif |
|---|
| 173 | totmasstmp=totmasstmp+masstmp(i) ! mass totale to be transferred if assume only one adjecent point |
|---|
| 174 | ENDDO |
|---|
| 175 | if (totmasstmp.gt.0.) then |
|---|
| 176 | !print*, 'totmass=',totmasstmp |
|---|
| 177 | !! 2) Available mass to be transferred |
|---|
| 178 | stock=maxval(masstmp(:))/secu ! kg/m2/s |
|---|
| 179 | !! 3) Real amounts of ice to be transferred : |
|---|
| 180 | masstmp(:)=masstmp(:)/totmasstmp*stock*area(ig)*iim/area(inddeb) !kg/m2/s all area around the pole are the same |
|---|
| 181 | !! 4) Tendencies |
|---|
| 182 | zdqsurf(ig)=-stock !kg/m2/s |
|---|
| 183 | !! 5) Update PQSURF |
|---|
| 184 | |
|---|
| 185 | ! move CH4 and CO too |
|---|
| 186 | if (methane) then |
|---|
| 187 | !! Variation of CH4 ice |
|---|
| 188 | 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) |
|---|
| 189 | !! remove CH4 |
|---|
| 190 | pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 |
|---|
| 191 | !! add CH4 |
|---|
| 192 | do i=inddeb,indfin |
|---|
| 193 | pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 |
|---|
| 194 | enddo |
|---|
| 195 | endif |
|---|
| 196 | if (carbox) then |
|---|
| 197 | !! Variation of CO ice |
|---|
| 198 | 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) |
|---|
| 199 | !! remove CO |
|---|
| 200 | pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco |
|---|
| 201 | !! add CO |
|---|
| 202 | do i=inddeb,indfin |
|---|
| 203 | pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco |
|---|
| 204 | enddo |
|---|
| 205 | endif |
|---|
| 206 | |
|---|
| 207 | ! Add N2 |
|---|
| 208 | do i=inddeb,indfin |
|---|
| 209 | pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep |
|---|
| 210 | !print*, 'height i=',phisfi0(i)/g+pqsurf(i,igcm_n2)/1000. |
|---|
| 211 | enddo |
|---|
| 212 | ! Remove N2 |
|---|
| 213 | pqsurf(ig,igcm_n2)=pqsurf(ig,igcm_n2)+zdqsurf(ig)*timstep |
|---|
| 214 | !print*, 'height ig=',phisfi0(ig)/g+pqsurf(ig,igcm_n2)/1000. |
|---|
| 215 | !print*, 'CHECK ig=',ig,zdqsurf(ig),sum(masstmp) |
|---|
| 216 | |
|---|
| 217 | endif |
|---|
| 218 | !print*, 'APPEL2',sum(masstmp),stock,totmasstmp,zdqsurf(ig),timstep,zdqsurf(ig)*timstep,sum(masstmp)*timstep |
|---|
| 219 | !call testconservmass2d(ngrid,pqsurf(:,1)) |
|---|
| 220 | |
|---|
| 221 | !!!!!!!!!!!!!!!!!!!! NOT THE POLE !!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 222 | ! here: ig = grid point with large amount of ice |
|---|
| 223 | ! Loop on each adjacent point |
|---|
| 224 | ! looking for adjacent point |
|---|
| 225 | ELSE |
|---|
| 226 | |
|---|
| 227 | !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point |
|---|
| 228 | DO compt=1,4 |
|---|
| 229 | i=edge(compt) |
|---|
| 230 | diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! (m) |
|---|
| 231 | if (diff.gt.difflim) then |
|---|
| 232 | !print*, 'diff=',diff,ig,i,compt,cas |
|---|
| 233 | masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) |
|---|
| 234 | else |
|---|
| 235 | masstmp(i)=0. |
|---|
| 236 | endif |
|---|
| 237 | totmasstmp=totmasstmp+masstmp(i) |
|---|
| 238 | ENDDO |
|---|
| 239 | if (totmasstmp.gt.0) then |
|---|
| 240 | !! 2) Available mass to be transferred |
|---|
| 241 | stock=maxval(masstmp(:))/secu ! kg/m2/s |
|---|
| 242 | !! 3) Real amounts of ice to be transferred : |
|---|
| 243 | DO compt=1,4 |
|---|
| 244 | i=edge(compt) |
|---|
| 245 | masstmp(i)=masstmp(i)/totmasstmp*stock*area(ig)/area(i) !kg/m2/s |
|---|
| 246 | !print*, 'masstmp=',masstmp(i) |
|---|
| 247 | ENDDO |
|---|
| 248 | !! 4) Tendencies |
|---|
| 249 | zdqsurf(ig)=-stock |
|---|
| 250 | !! 5) Update PQSURF |
|---|
| 251 | |
|---|
| 252 | ! move CH4 and CO too |
|---|
| 253 | if (methane) then |
|---|
| 254 | !! Variation of CH4 ice |
|---|
| 255 | 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) |
|---|
| 256 | !! remove CH4 |
|---|
| 257 | pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 |
|---|
| 258 | !! add CH4 |
|---|
| 259 | DO compt=1,4 |
|---|
| 260 | i=edge(compt) |
|---|
| 261 | pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 |
|---|
| 262 | enddo |
|---|
| 263 | endif |
|---|
| 264 | if (carbox) then |
|---|
| 265 | !! Variation of CO ice |
|---|
| 266 | 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) |
|---|
| 267 | !! remove CO |
|---|
| 268 | pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco |
|---|
| 269 | !! add CO |
|---|
| 270 | DO compt=1,4 |
|---|
| 271 | i=edge(compt) |
|---|
| 272 | pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco |
|---|
| 273 | enddo |
|---|
| 274 | endif |
|---|
| 275 | |
|---|
| 276 | ! Add N2 |
|---|
| 277 | DO compt=1,4 |
|---|
| 278 | i=edge(compt) |
|---|
| 279 | pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep |
|---|
| 280 | enddo |
|---|
| 281 | ! Remove N2 |
|---|
| 282 | pqsurf(ig,igcm_n2)=pqsurf(ig,igcm_n2)+zdqsurf(ig)*timstep |
|---|
| 283 | |
|---|
| 284 | endif |
|---|
| 285 | |
|---|
| 286 | if (pqsurf(ig,igcm_n2).lt.0) then |
|---|
| 287 | print*, pqsurf(ig,igcm_n2) |
|---|
| 288 | write(*,*) 'bug in spreadglacier' |
|---|
| 289 | stop |
|---|
| 290 | endif |
|---|
| 291 | |
|---|
| 292 | ENDIF ! cas |
|---|
| 293 | ENDIF ! qlim |
|---|
| 294 | |
|---|
| 295 | ENDDO ! ig |
|---|
| 296 | |
|---|
| 297 | end subroutine spreadglacier_paleo |
|---|
| 298 | |
|---|
| 299 | real function dist_pluto(lat1,lon1,lat2,lon2) |
|---|
| 300 | implicit none |
|---|
| 301 | #include "comcstfi.h" |
|---|
| 302 | real, intent(in) :: lat1,lon1,lat2,lon2 |
|---|
| 303 | real dlon,dlat |
|---|
| 304 | real a,c |
|---|
| 305 | real radi |
|---|
| 306 | radi=rad/1000. |
|---|
| 307 | |
|---|
| 308 | dlon = lon2 - lon1 |
|---|
| 309 | dlat = lat2 - lat1 |
|---|
| 310 | a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 |
|---|
| 311 | c = 2 * atan2( sqrt(a), sqrt(1-a) ) |
|---|
| 312 | dist_pluto = radi * c !! km |
|---|
| 313 | return |
|---|
| 314 | end function dist_pluto |
|---|
| 315 | |
|---|
| 316 | real function massflow(ig1,ig2,pts,dif,pqs,dt,slid) |
|---|
| 317 | !! Mass of ice to be transferred from one grid point to another |
|---|
| 318 | use comgeomfi_h |
|---|
| 319 | implicit none |
|---|
| 320 | #include "dimensions.h" |
|---|
| 321 | #include "dimphys.h" |
|---|
| 322 | #include "comcstfi.h" |
|---|
| 323 | #include "surfdat.h" |
|---|
| 324 | #include "comvert.h" |
|---|
| 325 | #include "callkeys.h" |
|---|
| 326 | #include "tracer.h" |
|---|
| 327 | real, intent(in) :: pts,dif,pqs,dt |
|---|
| 328 | integer, intent(in) :: ig1,ig2,slid |
|---|
| 329 | real lref,ac,n,theta,hdelta,ha,tau,qdry,qsl,usl0 |
|---|
| 330 | REAL dist_pluto ! function |
|---|
| 331 | |
|---|
| 332 | lref=dist_pluto(lati(ig2),long(ig2),lati(ig1),long(ig1))*1000. ! m |
|---|
| 333 | usl0=6.3e-6 ! m/s |
|---|
| 334 | ac=0.005*exp(9.37-422./pts) !0.005*exp(422./45.-422./pts) |
|---|
| 335 | n=2.1+0.0155*(pts-45.) |
|---|
| 336 | theta=atan(dif/(lref)) |
|---|
| 337 | hdelta=pts/20. |
|---|
| 338 | ha=pts*pts/8440. !pts*pts/20./422. |
|---|
| 339 | 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.) |
|---|
| 340 | qdry=qdry*0.5*exp( (pqs/1000./ha) / ( 1+ pqs /hdelta) ) |
|---|
| 341 | qsl=(pqs/1000.)/abs(tan(theta))*usl0*slid |
|---|
| 342 | tau=pqs/1000.*lref/abs((qdry+qsl)*tan(theta)) |
|---|
| 343 | massflow=pqs*(1.-exp(-dt/tau))/dt |
|---|
| 344 | return |
|---|
| 345 | end function massflow |
|---|
| 346 | |
|---|
| 347 | |
|---|
| 348 | |
|---|
| 349 | |
|---|
