subroutine spreadglacier_paleo(ngrid,nq,pqsurf, & phisfi0,timstep,ptsurf) use callkeys_mod, only: carbox, methane use comcstfi_mod, only: g use comgeomfi_h use geometry_mod, only: latitude, longitude, cell_area use radinc_h, only : naerkind use tracer_h, only: igcm_n2,igcm_ch4_ice,igcm_co_ice implicit none !================================================================== ! Purpose ! ------- ! Spreading the glacier : N2, cf Umurhan et al 2017 ! ! Inputs ! ------ ! ngrid Number of vertical columns ! pqsurf(ngrid,nq) N2 ice tracer on surface kg/m2 ! phisfi0 = phisfinew the geopotential of the bedrock ! below the N2 ice ! ptsurf surface temperature ! timstep dstep = timestep in pluto day for the call of glacial flow ! Outputs ! ------- ! pqsurf(ngrid,nq) new value for N2 ice tracer on surface kg/m2 ! ! Authors ! ------- ! Tanguy Bertrand (2015,2022) ! !================================================================== #include "dimensions.h" !----------------------------------------------------------------------- ! Arguments INTEGER ngrid, nq, ig, i REAL :: pqsurf(ngrid,nq) REAL :: phisfi0(ngrid) REAL :: ptsurf(ngrid) REAL timstep !----------------------------------------------------------------------- ! Local variables REAL tgla(ngrid),tbase(ngrid) !temperature at the base of glacier different of ptsurf REAL :: zdqsurf(ngrid) ! tendency of qsurf REAL massflow ! function REAL hlim,qlim,difflim,diff,stock,H0,totmasstmp INTEGER compt !compteur INTEGER slid ! option slid 0 or 1 INTEGER :: edge(4) ! index of the adjecent points, N,S,E,W INTEGER cas,inddeb,indfin REAL secu,dqch4,dqco REAL :: masstmp(ngrid) REAL :: masstot !---------------- INPUT ------------------------------------------------ difflim=0.5 ! limit height (m) between two adjacent point to start spreading the glacier zdqsurf(:)=0. !--------------- Dimensions ------------------------------------- ! distance between two adjacent latitude !distlim_lat=pi*rad/jjm/1000. ! km !distlim_diag=(distlim_lat*distlim_lat+distlim_lon*distlim_lon)**0.5 !! Threshold for ice thickness for which we activate the glacial flow hlim=10. !m qlim=hlim*1000. ! kg !! Security for not depleted all ice too fast in one timestep secu=4 !************************************* ! Loop over all points !************************************* DO ig=1,ngrid !if (ig.eq.ngrid) then ! print*, 'qpole=',pqsurf(ig,igcm_n2),qlim !endif !************************************* ! if significant amount of ice, where qsurf > qlim !************************************* IF (pqsurf(ig,igcm_n2).gt.qlim) THEN !************************************* ! temp glacier with gradient 15 K / km !************************************* ! surface temperature pts tgla(ig)=ptsurf(ig) ! temperature at the base (we neglect convection) tbase(ig)=tgla(ig)+20.*pqsurf(ig,igcm_n2)/1.e6 ! Umurhan et al. if (tbase(ig).gt.63.) then slid=1 ! activate slide else slid=0 endif !************************************* ! Selection of the adjacent index depending on the grid point !************************************* !! poles IF (ig.eq.1) then cas=0 inddeb=2 ! First adj grid point indfin=iim+1 ! Last adj grid point ELSEIF (ig.eq.ngrid) then cas=10 inddeb=ngrid-iim indfin=ngrid-1 !! 9 other cases: edges East, west, or in the center , at edges north south or in the center ELSEIF (mod(ig,iim).eq.2) then ! Edge east : N,S,W,E IF (ig.eq.2) then cas=1 edge = (/1, ig+iim,ig-1+iim,ig+1 /) ELSEIF (ig.eq.ngrid-iim) then cas=3 edge = (/ig-iim, ngrid,ig-1+iim,ig+1 /) ELSE cas=2 edge = (/ig-iim, ig+iim,ig-1+iim,ig+1 /) ENDIF ELSEIF (mod(ig,iim).eq.1) then ! Edge west IF (ig.eq.iim+1) then cas=7 edge = (/1,ig+iim,ig-1,ig+1-iim /) ELSEIF (ig.eq.ngrid-1) then cas=9 edge = (/ig-iim,ngrid,ig-1,ig+1-iim /) ELSE cas=8 edge = (/ig-iim,ig+iim,ig-1,ig+1-iim /) ENDIF ELSE IF ((ig.gt.2).and.(ig.lt.iim+1)) then cas=4 edge = (/1,ig+iim,ig-1,ig+1 /) ELSEIF ((ig.gt.ngrid-iim).and.(ig.lt.ngrid-1)) then cas=6 edge = (/ig-iim,ngrid,ig-1,ig+1 /) ELSE cas=5 edge = (/ig-iim,ig+iim,ig-1,ig+1 /) ENDIF ENDIF masstmp(:)=0. ! mass to be transferred totmasstmp=0. ! total mass to be transferred H0=phisfi0(ig)/g+pqsurf(ig,igcm_n2)/1000. ! height of ice at ig (m) !************************************* !!!!!!!!!!!!!!!!! POLE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !************************************* IF ((cas.eq.0).or.(cas.eq.10)) then ! Mean over all latitudes near the pole !call testconservmass2d(ngrid,pqsurf(:,1)) !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point DO i=inddeb,indfin diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! Height difference between ig and adjacent points (m) if (diff.gt.difflim) then masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) ! Mass to be transferred (kg/m2/s) else masstmp(i)=0. endif totmasstmp=totmasstmp+masstmp(i) ! mass totale to be transferred if assume only one adjecent point ENDDO if (totmasstmp.gt.0.) then !! 2) Available mass to be transferred stock=maxval(masstmp(:))/secu ! kg/m2/s !! 3) Real amounts of ice to be transferred : masstmp(:)=masstmp(:)/totmasstmp*stock*cell_area(ig)/cell_area(inddeb) !kg/m2/s all area around the pole are the same !! 4) Tendencies zdqsurf(ig)=-stock !kg/m2/s !! 5) Update PQSURF ! move CH4 and CO too if (methane) then !! Variation of CH4 ice 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) !! remove CH4 pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 !! add CH4 do i=inddeb,indfin pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 enddo endif if (carbox) then !! Variation of CO ice 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) !! remove CO pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco !! add CO do i=inddeb,indfin pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco enddo endif ! Add N2 do i=inddeb,indfin pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep enddo ! Remove N2 pqsurf(ig,igcm_n2)=pqsurf(ig,igcm_n2)+zdqsurf(ig)*timstep endif !call testconservmass2d(ngrid,pqsurf(:,1)) !!!!!!!!!!!!!!!!!!!! NOT THE POLE !!!!!!!!!!!!!!!!!!!!!!!! ! here: ig = grid point with large amount of ice ! Loop on each adjacent point ! looking for adjacent point ELSE !! 1) Compute Diff H0-H1 and mass to be transferred for each adjacent point DO compt=1,4 i=edge(compt) diff=H0-(phisfi0(i)/g+pqsurf(i,igcm_n2)/1000.) ! (m) if (diff.gt.difflim) then masstmp(i)=massflow(ig,i,tgla(ig),diff,pqsurf(ig,igcm_n2),timstep,slid) else masstmp(i)=0. endif totmasstmp=totmasstmp+masstmp(i) ENDDO if (totmasstmp.gt.0) then !! 2) Available mass to be transferred stock=maxval(masstmp(:))/secu ! kg/m2/s !! 3) Real amounts of ice to be transferred : DO compt=1,4 i=edge(compt) masstmp(i)=masstmp(i)/totmasstmp*stock*cell_area(ig)/cell_area(i) !kg/m2/s ENDDO !! 4) Tendencies zdqsurf(ig)=-stock !! 5) Update PQSURF ! move CH4 and CO too if (methane) then !! Variation of CH4 ice 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) !! remove CH4 pqsurf(ig,igcm_ch4_ice)=pqsurf(ig,igcm_ch4_ice)+dqch4 !! add CH4 DO compt=1,4 i=edge(compt) pqsurf(i,igcm_ch4_ice)=pqsurf(i,igcm_ch4_ice)-masstmp(i)/stock*dqch4 enddo endif if (carbox) then !! Variation of CO ice 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) !! remove CO pqsurf(ig,igcm_co_ice)=pqsurf(ig,igcm_co_ice)+dqco !! add CO DO compt=1,4 i=edge(compt) pqsurf(i,igcm_co_ice)=pqsurf(i,igcm_co_ice)-masstmp(i)/stock*dqco enddo endif ! Add N2 totmasstmp=0. DO compt=1,4 i=edge(compt) pqsurf(i,igcm_n2)=pqsurf(i,igcm_n2)+masstmp(i)*timstep totmasstmp=totmasstmp+masstmp(i)*cell_area(i) !! TB22 tmp enddo ! Remove N2 pqsurf(ig,igcm_n2)=pqsurf(ig,igcm_n2)+zdqsurf(ig)*timstep endif if (pqsurf(ig,igcm_n2).lt.0) then print*, pqsurf(ig,igcm_n2) write(*,*) 'bug in spreadglacier' stop endif ENDIF ! cas ENDIF ! qlim ENDDO ! ig end subroutine spreadglacier_paleo real function dist_pluto(lat1,lon1,lat2,lon2) use comcstfi_mod, only: rad implicit none real, intent(in) :: lat1,lon1,lat2,lon2 real dlon,dlat real a,c real radi radi=rad/1000. dlon = lon2 - lon1 dlat = lat2 - lat1 a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2 c = 2 * atan2( sqrt(a), sqrt(1-a) ) dist_pluto = radi * c !! km return end function dist_pluto real function massflow(ig1,ig2,pts,dif,pqs,dt,slid) !! Mass of ice to be transferred from one grid point to another use comgeomfi_h use geometry_mod, only: latitude, longitude, cell_area use comcstfi_mod, only: g, rad use tracer_h, only: rho_n2 implicit none #include "dimensions.h" real, intent(in) :: pts,dif,pqs,dt integer, intent(in) :: ig1,ig2,slid real lref,ac,n,theta,hdelta,ha,tau,qdry,qsl,usl0 REAL dist_pluto ! function lref=dist_pluto(latitude(ig2),longitude(ig2),latitude(ig1),longitude(ig1))*1000. ! m usl0=6.3e-6 ! m/s ac=0.005*exp(9.37-422./pts) !0.005*exp(422./45.-422./pts) n=2.1+0.0155*(pts-45.) theta=atan(dif/(lref)) hdelta=pts/20. ha=pts*pts/8440. !pts*pts/20./422. 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.) qdry=qdry*0.5*exp( (pqs/1000./ha) / ( 1+ pqs /hdelta) ) qsl=(pqs/1000.)/abs(tan(theta))*usl0*slid tau=pqs/1000.*lref/abs((qdry+qsl)*tan(theta)) massflow=pqs*(1.-exp(-dt/tau))/dt return end function massflow