Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

spreadglacier_paleo.F90_saved @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 18 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 13.8 KB

Rev	Line
[3175]	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_latdistlim_lat+distlim_londistlim_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(:)/totmasstmpstockarea(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)/totmasstmpstockarea(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.005exp(9.37-422./pts) !0.005exp(422./45.-422./pts)
	335	n=2.1+0.0155*(pts-45.)
	336	theta=atan(dif/(lref))
	337	hdelta=pts/20.
	338	ha=ptspts/8440. !ptspts/20./422.
	339	qdry=ac(rho_n2g1.e-6)n(pqs/1000.)*(n+2)/(n+2)tan(theta)*(n-1)/(1+tan(theta)tan(theta))**(n/2.)
	340	qdry=qdry0.5exp( (pqs/1000./ha) / ( 1+ pqs /hdelta) )
	341	qsl=(pqs/1000.)/abs(tan(theta))usl0slid
	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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format