Context Navigation

source: branches/iLoveclim/SOURCES/furst_schoof_mod.f90 @ 187

Last change on this file since 187 was 187, checked in by aquiquet, 6 years ago
Grisli-iloveclim branch merged to trunk at revision 185
File size: 26.3 KB

Rev	Line
[93]	1	! Module Furst Schoof
	2
	3	module furst_schoof_mod
	4
	5	use module3d_phy
	6	use deformation_mod_2lois
	7
	8	implicit none
	9
[110]	10	real :: frot_coef ! afq: fixed in Tsai but = beta in Schoof (for m_weert = 1)
[93]	11	real :: alpha_flot
	12	integer :: gr_select
[104]	13	real,dimension(nx,ny) :: back_force_x
	14	real,dimension(nx,ny) :: back_force_y
[110]	15	real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model???
[187]	16	real, parameter :: inv_mweert = 1./m_weert
[93]	17
	18	contains
	19
	20	subroutine init_furst_schoof
	21
	22	namelist/furst_schoof/frot_coef,gr_select
	23
	24	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
	25	read(num_param,furst_schoof)
	26
	27	write(num_rep_42,*)'!___________________________________________________________'
	28	write(num_rep_42,*)'! module furst_schoof'
	29	write(num_rep_42,*)
	30	write(num_rep_42,*)'! frot_coef = ',frot_coef
	31	write(num_rep_42,*)'! gr_select = ',gr_select
	32	write(num_rep_42,*)'! frot_coef : solid friction law coef 0.035 in Ritz et al 2001'
	33	write(num_rep_42,*)'! gr_select = 1 : Tsai , 2 : Schoof'
	34	write(num_rep_42,*)'!___________________________________________________________'
	35
	36	! frot = 0.035 ! f in the solid friction law, 0.035 in Ritz et al. 2001
[110]	37
[93]	38	! choix Tsai (loi de Coulomb) ou Schoof (loi de Weertman)
	39	alpha_flot= ro/row
[104]	40	back_force_x(:,:)=0.1
	41	back_force_y(:,:)=0.1
[93]	42
	43	end subroutine init_furst_schoof
	44
	45	subroutine interpol_glflux
	46
	47	implicit none
	48	real :: xpos,ypos
	49	real :: Hgl
	50	! integer :: igl
	51
	52	real :: phi_im2, phi_im1, phi_i, phi_ip1, phi_ip2, phi_ip3
	53	real :: phi_jm2, phi_jm1, phi_j, phi_jp1, phi_jp2, phi_jp3
	54	real,dimension(nx,ny) :: phi_xgl, phi_ygl ! flux gr line sur maille stag
[94]	55	integer,dimension(nx,ny) :: countx, county ! how often do we modify ux/uy
[93]	56	real :: phi_prescr
[94]	57	real :: toutpetit = 1e-6
	58	real :: denom, prodscal
[93]	59
[187]	60	real :: bfx, bfy
	61
[93]	62	!debug
	63	real,dimension(nx,ny) :: xpos_tab=9999.
	64	real,dimension(nx,ny) :: ypos_tab=9999.
	65	real,dimension(nx,ny) :: Hglx_tab=9999.
	66	real,dimension(nx,ny) :: Hgly_tab=9999.
	67	real,dimension(nx,ny) :: Hmxloc, Hmyloc ! pour eviter division par 0
	68	real,dimension(nx,ny) :: phi_prescr_tabx, phi_prescr_taby
	69
	70	real,dimension(nx,ny) :: uxbartemp, uybartemp
[96]	71
	72	real,dimension(nx,ny) :: archimtab
[104]	73
	74	real,dimension(nx,ny) :: imx_diag_in
[93]	75
	76	phi_prescr_tabx=0.
	77	phi_prescr_taby=0.
	78
	79	Hmxloc(:,:)= max(Hmx(:,:),1.)
	80	Hmyloc(:,:)=max(Hmy(:,:),1.)
	81
	82	uxbartemp(:,:)=0.
	83	uybartemp(:,:)=0.
[94]	84	countx(:,:)=0.
	85	county(:,:)=0.
[109]	86	gr_line_schoof(:,:)=0
[110]	87
[97]	88	archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel
[96]	89
[104]	90	imx_diag_in(:,:) = imx_diag(:,:)
	91
[93]	92	! detection des noeuds grounding line
	93	do j= 3,ny-3
	94	do i=3,nx-3
[97]	95	!archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
	96	if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded and with ice
[93]	97	! print*,'schoof gr line ij',i,j
	98	! selon x (indice i)
[96]	99
[110]	100	if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) then
[96]	101	!if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then ! grounding line between i-1,i CAS WEST (ecoulement vers grille West)
[93]	102
	103	call calc_xgl4schoof(-dx,alpha_flot,H(i,j),Bsoc(i,j),H(i-1,j),Bsoc(i-1,j),Sealevel,xpos,Hgl)
	104	xpos_tab(i,j)=xpos
	105	Hglx_tab(i,j)=Hgl
[187]	106
	107	! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
	108	if (xpos .lt. -dx/2.) then
	109	bfx = back_force_x(i,j)
	110	else
	111	bfx = back_force_x(i+1,j)
	112	endif
[93]	113	if (gr_select.eq.1) then ! flux de Tsai
[187]	114	call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
[110]	115	else if (gr_select.eq.2) then ! flux de Schoof
	116	! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
	117	if (xpos .lt. -dx/2.) then
	118	frot_coef = betamx(i,j)
	119	else
	120	frot_coef = betamx(i+1,j)
	121	endif
[187]	122	call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
[93]	123	else
[110]	124	print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
[93]	125	stop
	126	endif
	127	phi_prescr=-phi_prescr ! Doit Ãªtre negatif dans le cas West
	128	phi_prescr_tabx(i,j)=phi_prescr
	129
	130	if (xpos .lt. -dx/2.) then ! GL a west du i staggered, o centre, x stag
	131
	132	! grille ! .....x......o......x......o...G..x......O......x......o......x......o
	133	! ! i-2 i-1 i i+1 i+2
	134	!flux ! in out G out in
	135
	136
[110]	137	if (.not.(ilemx(i,j)).and..not.(ilemx(i-1,j))) then
	138
[93]	139	! extrapolation pour avoir uxbar(i-1,j) qui sera ensuite prescrit dans call rempli_L2
	140
	141	phi_im2 = Uxbar(i-2,j) * Hmxloc(i-2,j)
	142	phi_xgl(i-1,j) = phi_im2 + dx * (phi_prescr-phi_im2)/(2.5 * dx + xpos)
[94]	143	uxbartemp(i-1,j) = uxbartemp(i-1,j)+ phi_xgl(i-1,j) / Hmxloc(i-1,j) ! attention division par 0 possible ?
	144	countx (i-1,j) = countx(i-1,j)+1
	145	imx_diag (i-1,j) = 1
[109]	146	gr_line_schoof(i-1,j) = 1
[110]	147
[93]	148	! extrapolation pour avoir uxbar(i,j)
	149
	150	phi_ip1 = Uxbar(i+1,j) * Hmxloc(i+1,j)
	151	phi_xgl(i,j) = phi_ip1 + dx * (phi_prescr - phi_ip1)/(0.5 * dx - xpos)
	152	! print*,'uxbar(i,j)=',uxbar(i,j)
[94]	153	uxbartemp(i,j) = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j)
	154	countx (i,j) = countx(i,j)+1
	155	imx_diag (i,j) = 1
[109]	156	gr_line_schoof(i,j) = 1
[93]	157	! print*,'schoof uxbar(i,j)=',uxbar(i,j)
	158	! print*,'Hgl',Hgl
	159	! print*,'phi_prescr',phi_prescr, phi_prescr/Hgl
	160	! print*
	161
[110]	162	end if
	163
[93]	164	else ! GL a Est du i staggered, o centre, x stag
	165
	166	! grille ! .....x......o......x......o......x...G..O......x......o......x......o
	167	! i-2 i-1 i i+1 i+2
	168	!flux ! in in out G out in
	169
[110]	170	if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then
[93]	171
	172	! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
	173
	174	phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j)
	175	phi_xgl(i,j) = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos)
[94]	176	uxbartemp(i,j) = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ?
	177	countx (i,j) = countx(i,j)+1
	178	imx_diag (i,j) = 1
[110]	179	gr_line_schoof(i,j) = 1
	180
[93]	181	! extrapolation pour avoir uxbar(i+1,j)
	182
	183	phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j)
	184	phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5 * dx - xpos)
[94]	185	uxbartemp(i+1,j) = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j)
	186	countx (i+1,j) = countx(i+1,j)+1
	187	imx_diag (i+1,j) = 1
[109]	188	gr_line_schoof(i+1,j) = 1
[110]	189
	190	end if
	191
[93]	192	end if
[96]	193
[110]	194	else if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) then
[96]	195	!else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then ! grounding line between i,i+1 CAS EST (ecoulement vers grille Est)
[93]	196
	197	call calc_xgl4schoof(dx,alpha_flot,H(i,j),Bsoc(i,j),H(i+1,j),Bsoc(i+1,j),Sealevel,xpos,Hgl)
	198	xpos_tab(i,j)=xpos
	199	Hglx_tab(i,j)=Hgl
[187]	200
	201	! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
	202	if (xpos .lt. dx/2.) then
	203	bfx = back_force_x(i,j)
	204	else
	205	bfx = back_force_x(i+1,j)
	206	endif
[93]	207	if (gr_select.eq.1) then ! flux de Tsai
[187]	208	call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
[110]	209	else if (gr_select.eq.2) then ! flux de Schoof
	210	! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
[187]	211	if (xpos .lt. dx/2.) then
[110]	212	frot_coef = betamx(i,j)
	213	else
	214	frot_coef = betamx(i+1,j)
	215	endif
[187]	216	call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
[93]	217	else
[110]	218	print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
[93]	219	stop
	220	endif
	221	phi_prescr_tabx(i,j)=phi_prescr
[187]	222	if (xpos .lt. dx/2.) then ! GL a west du i staggered, o centre, x stag
[93]	223
	224	! grille ! .....x......o......x......o......x......O..G...x......o......x......o
	225	! ! i-2 i-1 i i+1 i+2
	226	!flux ! in out G out in
	227
[110]	228	if (.not.(ilemx(i,j)).and..not.(ilemx(i+1,j))) then
	229
[93]	230	! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
	231
	232	phi_im1 = Uxbar(i-1,j) * Hmxloc(i-1,j)
	233	phi_xgl(i,j) = phi_im1 + dx * (phi_prescr-phi_im1)/(1.5 * dx + xpos)
[94]	234	uxbartemp(i,j) = uxbartemp(i,j) + phi_xgl(i,j) / Hmxloc(i,j) ! attention division par 0 possible ?
	235	countx (i,j) = countx(i,j)+1
	236	imx_diag (i,j) = 1
[110]	237	gr_line_schoof(i,j) = 1
	238
[93]	239	! extrapolation pour avoir uxbar(i+1,j)
	240
	241	phi_ip2 = Uxbar(i+2,j) * Hmxloc(i+2,j)
	242	phi_xgl(i+1,j) = phi_ip2 + dx * (phi_prescr - phi_ip2)/(1.5* dx - xpos)
[94]	243	uxbartemp(i+1,j) = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j)
	244	countx (i+1,j) = countx(i+1,j)+1
	245	imx_diag (i+1,j) = 1
[110]	246	gr_line_schoof(i+1,j) = 1
	247
	248	end if
	249
[93]	250	else ! GL a west du i staggered, o centre, x stag
	251
	252	! grille ! ......x......o......x......O......x...G..o......x......o......x.....o
	253	! ! i-1 i i+1 i+2 i+3
	254	!flux ! in out G out in
	255
[110]	256	if (.not.(ilemx(i+1,j)).and..not.(ilemx(i+2,j))) then
	257
[93]	258	! extrapolation pour avoir uxbar(i,j) qui sera ensuite prescrit dans call rempli_L2
	259
	260	phi_i = Uxbar(i,j) * Hmxloc(i,j)
	261	phi_xgl(i+1,j) = phi_i + dx * (phi_prescr-phi_i)/(0.5 * dx + xpos)
[94]	262	uxbartemp(i+1,j) = uxbartemp(i+1,j) + phi_xgl(i+1,j) / Hmxloc(i+1,j) ! attention division par 0 possible ?
	263	countx (i+1,j) = countx(i+1,j)+1
	264	imx_diag (i+1,j) = 1
[110]	265	gr_line_schoof(i+1,j) = 1
[93]	266
	267	! extrapolation pour avoir uxbar(i+1,j)
	268
	269	phi_ip3 = Uxbar(i+3,j) * Hmxloc(i+3,j)
	270	phi_xgl(i+2,j) = phi_ip3 + dx * (phi_prescr - phi_ip3)/(2.5 * dx - xpos)
[94]	271	uxbartemp(i+2,j) = uxbartemp(i+2,j) + phi_xgl(i+2,j) / Hmxloc(i+2,j)
	272	countx (i+2,j) = countx(i+2,j)+1
	273	imx_diag (i+2,j) = 1
[109]	274	gr_line_schoof(i+2,j) = 1
[110]	275
	276	end if
	277
[93]	278	end if
	279
	280	end if
	281
	282	!*******************************************************************************
	283	! selon y (indice j)
[110]	284	if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) then
[96]	285	!if (flot(i,j-1).and.Uybar(i,j).LT.0.) then ! grounding line between j-1,j CAS SUD (ecoulement vers grille SUD)
[93]	286
	287
	288	call calc_xgl4schoof(-dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j-1),Bsoc(i,j-1),Sealevel,ypos,Hgl)
	289	ypos_tab(i,j)=ypos
	290	Hgly_tab(i,j)=Hgl
[187]	291
	292	! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
	293	if (ypos .lt. -dy/2.) then
	294	bfy = back_force_y(i,j)
	295	else
	296	bfy = back_force_y(i,j+1)
	297	endif
[93]	298	if (gr_select.eq.1) then ! flux de Tsai
[187]	299	call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
[110]	300	else if (gr_select.eq.2) then ! flux de Schoof
	301	! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
	302	if (ypos .lt. -dy/2.) then
	303	frot_coef = betamy(i,j)
	304	else
	305	frot_coef = betamy(i,j+1)
	306	endif
[187]	307	call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
[93]	308	else
[110]	309	print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
[93]	310	stop
	311	endif
	312	phi_prescr=-phi_prescr ! Doit Ãªtre negatif dans le cas Sud
	313	phi_prescr_taby(i,j)=phi_prescr
	314	if (ypos .lt. -dy/2.) then ! GL au sud du j staggered, o centre, x stag
	315
	316	! grille ! .....x......o......x......o...G..x......O......x......o......x......o
	317	! ! j-2 j-1 j j+1 j+2
	318	!flux ! in out G out in
	319
[110]	320	if (.not.(ilemy(i,j)).and..not.(ilemy(i,j-1))) then
[93]	321
	322	! extrapolation pour avoir uybar(i,j-1) qui sera ensuite prescrit dans call rempli_L2
	323
	324	phi_jm2 = Uybar(i,j-2) * Hmyloc(i,j-2)
	325	phi_ygl(i,j-1) = phi_jm2 + dy * (phi_prescr-phi_jm2)/(2.5 * dy + ypos)
[94]	326	uybartemp(i,j-1) = uybartemp(i,j-1) + phi_ygl(i,j-1) / Hmyloc(i,j-1) ! attention division par 0 possible ?
	327	county (i,j-1) = county(i,j-1)+1
[97]	328	imy_diag (i,j-1) = 1
[109]	329	gr_line_schoof(i,j-1) = 1
[93]	330
	331	! extrapolation pour avoir uybar(i,j)
	332
	333	phi_jp1 = Uybar(i,j+1) * Hmyloc(i,j+1)
	334	phi_ygl(i,j) = phi_jp1 + dy * (phi_prescr - phi_jp1)/(0.5 * dy - ypos)
[94]	335	uybartemp(i,j) = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j)
	336	county (i,j) = county(i,j)+1
	337	imy_diag (i,j) = 1
[109]	338	gr_line_schoof(i,j) = 1
[93]	339
[110]	340
	341	end if
	342
[93]	343	else ! GL au nord du j staggered, o centre, x stag
	344
	345	! grille ! .....x......o......x......o......x...G..O......x......o......x......o
	346	! j-2 j-1 j j+1 j+2
	347	!flux ! in in out G out in
	348
[110]	349	if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then
	350
[93]	351	! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
	352
	353	phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1)
	354	phi_ygl(i,j) = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos)
[94]	355	uybartemp(i,j) = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ?
	356	county (i,j) = county(i,j)+1
[93]	357	imy_diag(i,j) = 1
[109]	358	gr_line_schoof(i,j) = 1
[93]	359
	360	! extrapolation pour avoir uybar(i,j+1)
	361
	362	phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2)
	363	phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5 * dy - ypos)
[94]	364	uybartemp(i,j+1) = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1)
	365	county (i,j+1) = county(i,j+1)+1
	366	imy_diag (i,j+1) = 1
[109]	367	gr_line_schoof(i,j+1) = 1
[110]	368
	369	end if
	370
[93]	371	end if
[96]	372
[110]	373	else if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) then
[96]	374	!else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then ! grounding line between j,j+1 CAS NORD (ecoulement vers grille Nord)
[93]	375
	376	call calc_xgl4schoof(dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j+1),Bsoc(i,j+1),Sealevel,ypos,Hgl)
	377	ypos_tab(i,j)=ypos
	378	Hgly_tab(i,j)=Hgl
[187]	379
	380	! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
	381	if (ypos .lt. dy/2.) then
	382	bfy = back_force_y(i,j)
	383	else
	384	bfy = back_force_y(i,j+1)
	385	endif
[93]	386	if (gr_select.eq.1) then ! flux de Tsai
[187]	387	call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
[110]	388	else if (gr_select.eq.2) then ! flux de Schoof
	389	! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
[187]	390	if (ypos .lt. dy/2.) then
[110]	391	frot_coef = betamy(i,j)
	392	else
	393	frot_coef = betamy(i,j+1)
	394	endif
[187]	395	call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
[93]	396	else
[110]	397	print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
[93]	398	stop
	399	endif
	400	phi_prescr_taby(i,j)=phi_prescr
[187]	401	if (ypos .lt. dy/2.) then ! GL au sud du j staggered, o centre, x stag
[93]	402
	403	! grille ! .....x......o......x......o......x......O..G...x......o......x......o
	404	! ! j-2 j-1 j j+1 j+2
	405	!flux ! in out G out in
	406
[110]	407	if (.not.(ilemy(i,j)).and..not.(ilemy(i,j+1))) then
[93]	408
	409	! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
	410
	411	phi_jm1 = Uybar(i,j-1) * Hmyloc(i,j-1)
	412	phi_ygl(i,j) = phi_jm1 + dy * (phi_prescr-phi_jm1)/(1.5 * dy + ypos)
[94]	413	uybartemp(i,j) = uybartemp(i,j) + phi_ygl(i,j) / Hmyloc(i,j) ! attention division par 0 possible ?
	414	county (i,j) = county(i,j)+1
	415	imy_diag (i,j) = 1
[109]	416	gr_line_schoof(i,j) = 1
[93]	417
	418	! extrapolation pour avoir uybar(i,j+1)
	419
	420	phi_jp2 = Uybar(i,j+2) * Hmyloc(i,j+2)
	421	phi_ygl(i,j+1) = phi_jp2 + dy * (phi_prescr - phi_jp2)/(1.5* dy - ypos)
[94]	422	uybartemp(i,j+1) = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1)
	423	county (i,j+1) = county(i,j+1)+1
	424	imy_diag (i,j+1) = 1
[109]	425	gr_line_schoof(i,j+1) = 1
[93]	426
[110]	427	end if
	428
[93]	429	else ! GL au nord du j staggered, o centre, x stag
	430
	431	! grille ! ......x......o......x......O......x...G..o......x......o......x.....o
	432	! ! j-1 j j+1 j+2 j+3
	433	!flux ! in out G out in
	434
[110]	435	if (.not.(ilemy(i,j+1)).and..not.(ilemy(i,j+2))) then
[93]	436
	437	! extrapolation pour avoir uybar(i,j) qui sera ensuite prescrit dans call rempli_L2
	438
	439	phi_j = Uybar(i,j) * Hmyloc(i,j)
	440	phi_ygl(i,j+1) = phi_j + dy * (phi_prescr-phi_j)/(0.5 * dy + ypos)
[94]	441	uybartemp(i,j+1) = uybartemp(i,j+1) + phi_ygl(i,j+1) / Hmyloc(i,j+1) ! attention division par 0 possible ?
	442	county (i,j+1) = county(i,j+1)+1
	443	imy_diag (i,j+1) = 1
[110]	444	gr_line_schoof(i,j+1) = 1
	445
[93]	446	! extrapolation pour avoir uybar(i,j+1)
	447
	448	phi_jp3 = Uybar(i,j+3) * Hmyloc(i,j+3)
	449	phi_ygl(i,j+2) = phi_jp3 + dy * (phi_prescr - phi_jp3)/(2.5 * dy - ypos)
[94]	450	uybartemp(i,j+2) = uybartemp(i,j+2) + phi_ygl(i,j+2) / Hmyloc(i,j+2)
	451	county (i,j+2) = county(i,j+2)+1
	452	imy_diag (i,j+2) = 1
[109]	453	gr_line_schoof(i,j+2) = 1
[110]	454
	455	end if
	456
[93]	457	end if
	458
	459	end if
	460	end if
	461	end do
	462	end do
[104]	463
	464	! afq -- if we discard the points with multiple fluxes coming, uncomment:
[110]	465	where (countx(:,:).ge.2)
	466	uxbartemp(:,:)=uxbar(i,j)*countx(:,:)
	467	imx_diag(:,:) = imx_diag_in(:,:)
	468	end where
	469	where (county(:,:).ge.2)
	470	uybartemp(:,:)=uybar(:,:)*county(:,:)
	471	imy_diag(:,:) = imx_diag_in(:,:)
	472	end where
[104]	473	! afq --
[94]	474
[104]	475	where (countx(:,:).ne.0) uxbar(:,:)= max(min(uxbartemp(:,:)/countx(:,:),1000.),-1000.)
	476	where (county(:,:).ne.0) uybar(:,:)= max(min(uybartemp(:,:)/county(:,:),1000.),-1000.)
	477
[94]	478	do j= 3,ny-3
	479	do i=3,nx-3
	480	denom = sqrt((uxbartemp(i,j)2 + uybartemp(i,j)2))sqrt((uxbar(i,j)2 + uybar(i,j)*2))
	481	if (denom.gt.toutpetit) then
	482	! prodscal is cos theta in : U . V = u v cos theta
	483	! this number is between -1 and 1, 1 = same direction & -1 = opposite direction
	484	! we could use this as an artificial back force...
	485	prodscal = (uxbartemp(i,j)uxbar(i,j) + uybartemp(i,j)uybar(i,j))/ &
	486	(sqrt((uxbartemp(i,j)2 + uybartemp(i,j)2))sqrt((uxbar(i,j)2 + uybar(i,j)*2)))
	487	endif
	488	end do
	489	end do
	490
[93]	491	!!$ do j=1,ny
	492	!!$ do i=1,nx
	493	!!$ write(598,*) xpos_tab(i,j)
	494	!!$ write(599,*) ypos_tab(i,j)
	495	!!$ write(600,*) Hglx_tab(i,j)
	496	!!$ write(601,*) Hgly_tab(i,j)
	497	!!$ write(602,*) Abar(i,j)
	498	!!$ write(603,*) phi_prescr_tabx(i,j)/Hglx_tab(i,j)
	499	!!$ write(604,*) phi_prescr_taby(i,j)/Hgly_tab(i,j)
	500	!!$ enddo
	501	!!$ enddo
[187]	502
	503	debug_3D(:,:,66) = phi_prescr_tabx(:,:)
	504	debug_3D(:,:,67) = phi_prescr_taby(:,:)
	505
[93]	506	end subroutine interpol_glflux
	507
	508	!-------------------------------------------------------------------------
	509	! calcule la position sous maille de la ligne dechouage
	510	subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,H_1,B_1,SL,xpos,Hgl)
	511
	512
	513	implicit none
	514
	515	real,intent(in) :: dy !< longueur orientee de la maille
	516	real,intent(in) :: alpha !< coefficient de flottaison = rho/rhow
	517	real,intent(in) :: H_0 !< epaisseur au point pose
	518	real,intent(in) :: B_0 !< altitude socle au point pose
	519	real,intent(in) :: H_1 !< epaisseur au point flottant
	520	real,intent(in) :: B_1 !< altitude socle au point flottant
	521	real,intent(in) :: SL !< Sea level
	522	real,intent(out) :: xpos !< position de la ligne (en distance depuis le point pose)
	523	real,intent(out) :: Hgl !< Epaisseur de glace a la grounding line
	524
	525
	526
	527	real :: Cgl ! variable de travail
	528
[96]	529	!~ Cgl = (-alpha * (H_1-H_0) - (B_1-B_0)) ! -(rho/rhow (H1-H0) + (B1-B0))
[93]	530
[96]	531	!~ if (abs(Cgl).gt.1.e-5) then
	532	!~ xpos = (B_0 + alpha * H_0 - SL) / Cgl
	533	!~ else
	534	!~ xpos = 1. - 1./abs(dy) ! verifier
	535	!~ end if
	536
	537	!~ if (xpos.LT.0..OR.xpos.GT.1.) then
	538	!~ print*,'calc_xgl4schoof : xpos value error i,j=',i,j
	539	!~ print*, 'xpos,Cgl', xpos,Cgl
	540	!~ print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
	541	!~ stop
	542	!~ endif
	543
	544	Cgl = (B_1-B_0) + alpha * (H_1-H_0)
	545
[93]	546	if (abs(Cgl).gt.1.e-5) then
[96]	547	xpos = (SL - (B_0 + alpha * H_0)) / Cgl
[93]	548	else
[96]	549	xpos = 0.5 ! verifier
[93]	550	end if
	551
[96]	552	! print*
	553	! print*, 'xpos', xpos
	554	! print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
	555	! print*,'B_1, H_1', B_1, H_1
	556
	557
[93]	558	if (xpos.LT.0..OR.xpos.GT.1.) then
[104]	559	print*,'calc_xgl4schoof : xpos value error i,j=',i,j
	560	print*, 'xpos,Cgl', xpos,Cgl
	561	print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
	562	print,'archim=',B_1+H_1alpha - SL
	563	!stop
	564	xpos = min(max(0.,xpos),1.)
[97]	565	endif
[96]	566
[93]	567	xpos = xpos * dy
	568
	569	Hgl= H_0 + (H_1 - H_0) * xpos/dy
	570
	571	end subroutine calc_xgl4schoof
	572
	573
	574	!--------------------------------------------------------------------
	575	! Schoof_flux : calcule le coefficient et le flux de Schoof
	576	!--------------------------------------------------------------------
	577
	578	! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement)
	579
	580	subroutine flux_Schoof4Schoof (Hgl,A_mean,C_frot,alpha,n_glen,m_weert,back_force_coef,phi_schoof)
	581
	582
	583	implicit none
	584
	585	real,intent(in) :: Hgl
	586	real,intent(in) :: A_mean
	587	real,intent(in) :: C_frot
	588	real,intent(in) :: alpha
	589	real,intent(in) :: n_glen
	590	real,intent(in) :: m_weert
	591	real,intent(in) :: back_force_coef
	592	real,intent(out) :: phi_schoof
	593
[110]	594	! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta
[93]	595
[187]	596	!phi_schoof = (((A_mean * rog*(n_glen+1) (1 - alpha)n_glen) / (4n_glen C_frot)) (1./(1.+1./m_weert)))back_force_coef**(n_glen/(1+1./m_weert))
	597	!phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert))
	598	phi_schoof = (((A_mean * rog*(n_glen+1.) (1. - alpha)n_glen) / (4.n_glenC_frot)) (1./(1.+inv_mweert)))back_force_coef**(n_glen/(1.+inv_mweert))
	599	phi_schoof = phi_schoof * Hgl**((n_glen+3.+inv_mweert)/(1.+inv_mweert))
[110]	600
[93]	601
	602
	603	end subroutine flux_Schoof4Schoof
	604
	605	!--------------------------------------------------------------------
	606	! Tsai_flux : calcule le coefficient et le flux de Schoof
	607	!--------------------------------------------------------------------
	608
	609	! Dans un premier temps, pas de terme horizontal (regarder dans toy_recul une tentative de confinement)
	610
	611	subroutine flux_Tsai4Schoof (Hgl,A_mean,f_frot,alpha,n_glen,back_force_coef,phi_Tsai)
	612
	613	implicit none
	614
	615	real,intent(in) :: Hgl
	616	real,intent(in) :: A_mean
	617	real,intent(in) :: f_frot
	618	real,intent(in) :: alpha
	619	real,intent(in) :: n_glen
	620	real,intent(in) :: back_force_coef
	621	real,intent(out) :: phi_Tsai
	622
	623	real, parameter :: Q0 = 0.61
	624
	625
[187]	626	phi_Tsai = Q0 * ((8. * A_mean * rog*(n_glen) (1. - alpha)(n_glen-1.)) / (4.n_glen f_frot))back_force_coef**(n_glen-1.)
	627	phi_Tsai = phi_Tsai * Hgl**(n_glen+2.)
[93]	628
	629	end subroutine flux_Tsai4Schoof
	630
	631	end module furst_schoof_mod

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

Download in other formats: