Context Navigation

source: trunk/SOURCES/flottab2-0.7.f90 @ 223

Last change on this file since 223 was 223, checked in by aquiquet, 6 years ago
Bug correction for S variation due to isostasy: S can rise back above sealevel after it has been depressed below sealevel
File size: 42.6 KB

Rev	Line
[4]	1	!> \file flottab2-0.7.f90
	2	!! Module pour determiner les endroits ou la glace
	3	!! flotte , les iles, et la position du front de glace
	4	!<
	5
	6	!> \namespace flottab_mod
	7	!! Determine les endroits ou la glace
	8	!! flotte , les iles, et la position du front de glace
	9	!! \author Vincent & Cat
	10	!! \date 10 juillet 2005
	11	!! @note Used module
	12	!! @note - use module3D_phy
	13	!! @note - use module_choix
	14	!! @todo "nasty Island". If the bedrock is above sealevel force grounded (mk_init)
	15	!<
	16	module flottab_mod
	17
[71]	18	!$ USE OMP_LIB
[4]	19	USE module3D_phy
	20	use module_choix
	21
	22
	23
	24	IMPLICIT NONE
	25
	26	real :: surnet !< surnet hauteur de glace au dessus de la mer
	27	real :: archim !< test de flottaison
[99]	28	! real, parameter :: Hmin=1.001 !< Hmin pour être considere comme point ice
[4]	29
	30	integer:: itestf
	31
	32	logical,dimension(nx,ny) :: gz1mx,gz1my
	33	logical,dimension(nx,ny) :: fl1mx,fl1my
	34
	35
	36	real,dimension(nx,ny) :: uxs1 !< uxbar a l'entree de flottab
	37	real,dimension(nx,ny) :: uys1 !< uybar a l'entree de flottab
	38
	39
	40	integer pmx,pmy !pm=plus-moins -1 ou 1 pour x et y
	41
	42
	43	! Variables pour la determination des differents shelfs/stream
	44	! (representés comme des taches ou l'on resoud l'eq elliptique)
	45	!________________________________________________________________
	46	integer,parameter :: n_ta_max=2000!< nombre de tache max
	47	integer,dimension(nx,ny) :: table_out !< pour les numeros des taches
	48	integer,dimension(nx,ny) :: tablebis !< pour les numeros des taches
	49	integer,dimension(0:n_ta_max) :: compt !< contient les equivalence entre les taches
	50	integer,dimension(0:n_ta_max) :: nb_pts_tache !< indique le nombre de points par tache
[111]	51	logical,dimension(0:n_ta_max) :: iceberg1D !< T si iceberg, F si calotte posee
[110]	52
[4]	53	logical,dimension(nx,ny) :: mask_tache_ij !< masque de travail
	54	!< vrai pour toute la tache de i,j
	55	integer,dimension(2) :: smax_coord !< pour le maxloc des iles
	56
	57	! Variables pour determiner le point le plus haut (surf)
	58	! d'une ile completement stream
	59	!_________________________________________________________
	60
	61	! icetrim : T si ice stream, F si calotte posee(vertical shear)
	62	logical,dimension(n_ta_max) :: icetrim
	63
	64	integer :: ii,jj
	65	integer :: smax_i
	66	integer :: smax_j
	67	real :: smax_
	68	integer :: numtache
	69	integer :: nb_pt
[142]	70	real :: petit_H=0.001 ! pour test ice sur zone flottante
[4]	71	contains
	72	! -----------------------------------------------------------------------------------
	73	!> SUBROUTINE: flottab()
	74	!! Cette routine determine les endroits ou la glace
	75	!! flotte , les iles, et la position du front de glace
	76	!! @note Il y a 4 sortes de zone
	77	!! @note - Pose
	78	! @note - Grounding zone et streams gzmx et gzmy
	79	! @note - Iles ilemx, ilemy
	80	! @note - flottant flot sur le noeud majeur, flotmx sur le noeud mineur
	81	!>
[86]	82	subroutine flottab()
	83	!
	84	! Vince 5 Jan 95
	85	! Modifie 20 Jan 95
	86	! Modifie le 30 Novembre 98
	87	! Passage f90 + determination des fronts Vincent dec 2003
	88	! nettoyage et nouvelle détermination de gzmx et gzmy Cat le 10 juillet 2005
	89	! Re-nettoyage par Cat en aout 2006.
	90	! Le calcul de gzmx et gzmy pour les points intérieurs passe dans la subroutine dragging
	91	! pour les points cotiers, toujours fait dans flottab
	92	!
	93	! -----------------
	94	!
	95	! Cette routine determine les endroits ou la glace
	96	! flotte , les iles, et la position du front de glace
	97	!
	98	! Il y a 4 sortes de zone
	99	! Pose
	100	! Grounding zone et streams gzmx et gzmy
	101	! Iles ilemx, ilemy
	102	! flottant flot sur le noeud majeur, flotmx sur le noeud mineur
[4]	103
[86]	104	! passage dans flottab tous les pas de temps dt !
	105	!
	106	! _________________________________________________________
[4]	107
[86]	108	!~ print*,'debut flottab',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel
	109	!~ print*,'debut flottab',flot(132,183),ice(132,183)
[99]	110	!print*,'H(90,179) flottab 1',H(90,179),ice(90,179), flot(90,179)
[4]	111
[86]	112	if (itracebug.eq.1) call tracebug(' Entree dans routine flottab')
[4]	113
[86]	114	SHELFY = .FALSE.
[4]	115
	116
[86]	117	! cas particulier des runs paleo ou on impose un masque grounded
[4]	118
[86]	119	!$OMP PARALLEL PRIVATE(archim,surnet)
	120	if (igrdline.eq.2) then
	121	!$OMP WORKSHARE
	122	where ( mk_init(:,:).eq.1) ! pose
	123	flot(:,:) = .False.
	124	H(:,:)=max(H(:,:),(10.+sealevel-Bsoc(:,:))*row/ro) ! pour avoir archim=10
	125	S(:,:) = Bsoc(:,:) + H(:,:)
	126	elsewhere
	127	flot(:,:) = .True.
	128	end where
	129	!$OMP END WORKSHARE
	130	end if
[4]	131
	132
[86]	133	! 1-INITIALISATION
	134	! ----------------
	135	! initialisation des variables pour detecter les points qui se mettent
	136	! a flotter entre 2 dtt
[4]	137
[86]	138	appel_new_flot=.false.
	139	!$OMP DO
	140	do j=1,ny
	141	do i=1,nx
	142	new_flot_point(i,j)=.false.
	143	new_flotmx(i,j)=.false.
	144	new_flotmy(i,j)=.false.
[4]	145	enddo
[86]	146	enddo
	147	!$OMP END DO
[4]	148
[86]	149	! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m
[4]	150
[86]	151	!$OMP WORKSHARE
	152	ICE(:,:)=0
	153	front(:,:)=0
	154	frontfacex(:,:)=0
	155	frontfacey(:,:)=0
[111]	156	isolx(:,:)=.false.
	157	isoly(:,:)=.false.
[86]	158	cotemx(:,:)=.false.
	159	cotemy(:,:)=.false.
	160	boost=.false.
[111]	161	iceberg(:,:)=.false.
[86]	162	!$OMP END WORKSHARE
[4]	163
[86]	164	! fin de l'initialisation
	165	!_____________________________________________________________________
[4]	166
[86]	167	! 2-TESTE LES NOUVEAUX POINTS FLOTTANTS
	168	! -------------------------------------
[4]	169
[86]	170	!$OMP DO
	171	do j=1,ny
	172	do i=1,nx
[4]	173
[86]	174	uxs1(i,j)=uxbar(i,j)
	175	uys1(i,j)=uybar(i,j)
[4]	176
[86]	177	archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
	178	! if ((i.eq.132).and.(j.eq.183)) print*,'archim=',archim
[4]	179
	180
[86]	181	arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.e-3)) then ! le point flotte
	182	mk(i,j)=1
[4]	183
	184
[86]	185	ex_pose: if ((.not.FLOT(I,J)).and.(isynchro.eq.1)) then ! il ne flottait pas avant
	186	FLOT(I,J)=.true.
	187	BOOST=.false.
[4]	188
[86]	189	if (igrdline.eq.1) then ! en cas de grounding line prescrite
	190	flot(i,j)=.false.
	191	H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro ! pour avoir archim=10
	192	new_flot_point(i,j)=.false.
	193	endif
[4]	194
[86]	195	else ! isynchro=0 ou il flottait déja
[4]	196
[86]	197	if (.not.FLOT(I,J)) then ! il ne flottait pas (isynchro=0)
	198	new_flot_point(i,j)=.true. ! signale un point qui ne flottait pas
	199	! au pas de temps precedent
	200	flot(i,j)=.true.
[4]	201
[86]	202	if (igrdline.eq.1) then ! en cas de grounding line prescrite
	203	flot(i,j)=.false.
	204	H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro ! pour avoir archim=10
	205	new_flot_point(i,j)=.false.
	206	endif
[4]	207
[86]	208	endif
	209	endif ex_pose
[4]	210
	211
[223]	212	else if ((H(i,j).ge.0.).and.(archim.GE.0.)) then ! le point ne flotte pas et est englace
[86]	213	mk(i,j)=0
[4]	214
[86]	215	if(FLOT(I,J)) then ! mais il flottait avant
[4]	216	FLOT(I,J)=.false.
	217	BOOST=.false.
[86]	218	endif
	219	!cdc correction topo pour suivre variations sealevel
	220	!cdd S(i,j)=Bsoc(i,j)+H(i,j)
	221	S(i,j)=Bsoc(i,j)+H(i,j)
	222	B(i,j)=Bsoc(i,j)
[4]	223
[86]	224	else if ((H(i,j).LE.0.).and.(archim.LT.0.)) then ! terre deglace qui devient ocean
[99]	225	!cdc ice(i,j)=0
	226	!cdc H(i,j)=1.
[102]	227	!cdc 1m H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01))
[121]	228	flot(i,j)=.true. !cdc points ocean sont flot meme sans glace
[102]	229	H(i,j)=0.
[119]	230	S(i,j)=sealevel
[86]	231	B(i,j)=S(i,j)-H(i,j)
	232	endif arch
	233
	234	! Si la glace flotte -> PRUDENCE !!!
	235	! S et B sont alors determines avec la condition de flottabilite
	236
	237	if (flot(i,j)) then
	238	shelfy = .true.
	239
	240	surnet=H(i,j)*(1.-ro/row)
	241	S(i,j)=surnet+sealevel
	242	B(i,j)=S(i,j)-H(i,j)
	243	end if
	244
	245	end do
[4]	246	end do
[86]	247	!$OMP END DO
	248	!~ print*,'flottab 2',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel
	249	!~ print*,'flottab 2',flot(132,183),ice(132,183)
	250	!~ if (S(132,183).LT.sealevel) print*,'BUGGGGGGGGGGGGG !!!!!!!!!!!!!!!!'
[119]	251
[4]	252	!!$ do i=1,nx
	253	!!$ do j=1,ny
	254	!!$ if (flot(i,j)) then
	255	!!$ mk(i,j)=1
	256	!!$ else
	257	!!$ mk(i,j)=0
	258	!!$ endif
	259	!!$ end do
	260	!!$ end do
	261
	262
	263
	264
[86]	265	!-----------------------------------------------------------------------
	266	!$OMP DO
	267	domain_x: do j=1,ny
	268	do i=2,nx
[4]	269
[86]	270	! 3_x A- NOUVELLE DEFINITION DE FLOTMX, LES POINTS
	271	! AYANT UN DES VOISINS FLOTTANTS SONT FLOTMX ET GZMX
	272	! -------------------------------------------
[4]	273
[86]	274	! fl1 est l'ancienne valeur de flotmx
	275	gz1mx(i,j)=gzmx(i,j) ! gz1 est l'ancienne valeur de gzmx
	276	fl1mx(i,j)=flotmx(i,j) ! fl1 est l'ancienne valeur de flotmx
[4]	277
[86]	278	flotmx(i,j)=flot(i,j).and.flot(i-1,j)
[4]	279
[86]	280	! test pour detecter les nouveaux flotmx entre 2 dtt :
[4]	281
[86]	282	if (flotmx(i,j).and.(new_flot_point(i,j).or. &
	283	new_flot_point(i-1,j))) then
	284	appel_new_flot=.true.
	285	new_flotmx(i,j)=.true.
	286	endif
[4]	287
	288
[86]	289	! premiere determination de gzmx
	290	!__________________________________________________________________________
[4]	291
[86]	292	! gzmx si un des deux voisins est flottant et l'autre posé
	293	! i-1 i
	294	gzmx(i,j)=((flot(i,j).and..not.flot(i-1,j)) & ! F P
	295	.or.(.not.flot(i,j).and.flot(i-1,j))) ! P F
[4]	296
[86]	297	! A condition d'etre assez proche de la flottaison
	298	! sur le demi noeud condition archim < 100 m
[4]	299
[86]	300	archim=(Bsoc(i,j)+Bsoc(i-1,j))0.5-sealevel+ro/rowHmx(i,j)
	301	gzmx(i,j)=gzmx(i,j).and.(archim.le.100.)
	302	cotemx(i,j)=gzmx(i,j)
[4]	303
[86]	304	end do
	305	end do domain_x
	306	!$OMP END DO
	307	!if (itracebug.eq.1) call tracebug(' routine flottab apres domain_x')
[4]	308
[86]	309	! 3_y B- NOUVELLE DEFINITION DE FLOTMY
	310	! --------------------------------
	311	!$OMP DO
	312	domain_y: do j=2,ny
	313	do i=1,nx
[4]	314
[86]	315	gz1my(i,j)=gzmy(i,j) ! gz1 est l'ancienne valeur de gzmy
	316	fl1my(i,j)=flotmy(i,j) ! fl1 est l'ancienne valeur de flotmy
[4]	317
[86]	318	flotmy(i,j)=flot(i,j).and.flot(i,j-1)
[4]	319
[86]	320	! test pour detecter les nouveaux flotmy entre 2 dtt :
[4]	321
[86]	322	if (flotmy(i,j).and.(new_flot_point(i,j).or. &
	323	new_flot_point(i,j-1))) then
	324	appel_new_flot=.true.
	325	new_flotmy(i,j)=.true.
	326	endif
[4]	327
[86]	328	! premiere determination de gzmy
	329	!__________________________________________________________________________
[4]	330
[86]	331	! gzmy si un des deux voisins est flottant et l'autre posé
[4]	332
[86]	333	gzmy(i,j)=((flot(i,j).and..not.flot(i,j-1)) &
	334	.or.(.not.flot(i,j).and.flot(i,j-1)))
[4]	335
[86]	336	! A condition d'etre assez proche de la flottaison
	337	! sur le demi noeud condition archim > 100 m
[4]	338
[86]	339	archim=(Bsoc(i,j)+Bsoc(i,j-1))0.5-sealevel+ro/rowHmy(i,j)
	340	gzmy(i,j)=gzmy(i,j).and.(archim.le.100.)
	341	cotemy(i,j)=gzmy(i,j)
[4]	342
[86]	343	end do
	344	end do domain_y
	345	!$OMP END DO
[4]	346
	347
[86]	348	!-------------------------------------------------------------------------------------
	349	! attention : pour expériences Heino
	350	! gzmy(i,j)=gzmy_heino(i,j)
	351	! shelfy=.true.
	352	! shelfy=.false.
	353	!____________________________________________________________________________________
[4]	354
	355
	356	!!$
	357	!!$
	358	!!$! 4- Condition sur les bords
	359	!!$
	360	!!$
	361	!!$ do i=2,nx
	362	!!$ flotmx(i,1) = (flot(i,1).or.flot(i-1,1))
	363	!!$ flotmy(i,1) = .false.
	364	!!$ end do
	365	!!$
	366	!!$ do j=2,ny
	367	!!$ flotmy(1,j) = (flot(1,j).or.flot(1,j-1))
	368	!!$ flotmx(1,j) = .false.
	369	!!$ end do
	370	!!$
	371	!!$ flotmx(1,1) = .false.
	372	!!$ flotmy(1,1) = .false.
	373	!!$
	374
	375
[86]	376	! 4- determination des iles
	377	! -------------------------
	378	!$OMP WORKSHARE
	379	ilemx(:,:)=.false.
	380	ilemy(:,:)=.false.
	381	!$OMP END WORKSHARE
[4]	382
[86]	383	! selon x
	384	!$OMP DO
	385	ilesx: do j=2,ny-1
	386	do i=3,nx-2
	387	! F G F
	388	! x
	389	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
	390	if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. &
	391	(sdx(i,j).LT.1.E-02)) then
[4]	392	ilemx(i,j)=.true.
	393	ilemx(i+1,j)=.true.
	394
[86]	395	! F G G F
	396	! x
	397	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
	398	else if ((flot(i-1,j).and..not.flot(i,j) &
	399	.and..not.flot(i+1,j)).and.flot(i+2,j).and. &
	400	(sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then
	401	ilemx(i,j)=.true.
	402	ilemx(i+1,j)=.true.
	403	ilemx(i+2,j)=.true.
[4]	404
[86]	405	! F G G F
	406	! x
	407	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
	408	else if ((flot(i-2,j).and..not.flot(i-1,j) &
	409	.and..not.flot(i,j)).and.flot(i+1,j).and. &
	410	(sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then
	411	ilemx(i-1,j)=.true.
	412	ilemx(i,j)=.true.
	413	ilemx(i+1,j)=.true.
[4]	414
[86]	415	! F G G G F
	416	! x
	417	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
	418	else if ((i.lt.nx-2) &
	419	.and.(flot(i-2,j).and..not.flot(i-1,j) &
	420	.and..not.flot(i,j)).and..not.flot(i+1,j) &
	421	.and.flot(i+2,j).and. &
	422	(sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02 &
	423	.and.sdx(i+1,j).LT.1.E-02)) then
	424	ilemx(i-1,j)=.true.
	425	ilemx(i,j)=.true.
	426	ilemx(i+1,j)=.true.
	427	ilemx(i+2,j)=.true.
[4]	428
[86]	429	endif
[4]	430
[86]	431	end do
	432	end do ilesx
	433	!$OMP END DO
[4]	434
[86]	435	! selon y
	436	!$OMP DO
	437	ilesy: do j=3,ny-2
	438	do i=2,nx-1
	439	! F G F
	440	! x
	441	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
	442	if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. &
	443	(sdy(i,j).LT.1.E-02)) then
[4]	444	ilemy(i,j)=.true.
	445	ilemy(i,j+1)=.true.
	446
[86]	447	! F G G F
	448	! x
	449	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
	450	else if ((flot(i,j-1).and..not.flot(i,j) &
	451	.and..not.flot(i,j+1)).and.flot(i,j+2).and. &
	452	(sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then
	453	ilemy(i,j)=.true.
	454	ilemy(i,j+1)=.true.
	455	ilemy(i,j+2)=.true.
[4]	456
[86]	457	! F G G F
	458	! x
	459	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
	460	else if ((flot(i,j-2).and..not.flot(i,j-1) &
	461	.and..not.flot(i,j)).and.flot(i,j+1).and. &
	462	(sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then
	463	ilemy(i,j-1)=.true.
	464	ilemy(i,j)=.true.
	465	ilemy(i,j+1)=.true.
[4]	466
[86]	467	! F G G G F
	468	! x
	469	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
	470	else if ((j.lt.ny-2) &
	471	.and.(flot(i,j-2).and..not.flot(i,j-1) &
	472	.and..not.flot(i,j)).and..not.flot(i,j+1) &
	473	.and.flot(i,j+2).and. &
	474	(sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02 &
	475	.and.sdy(i,j+1).LT.1.E-02)) then
	476	ilemy(i,j-1)=.true.
	477	ilemy(i,j)=.true.
	478	ilemy(i,j+1)=.true.
	479	ilemy(i,j+2)=.true.
	480	endif
	481	end do
	482	end do ilesy
[71]	483	!$OMP END DO
	484	!$OMP END PARALLEL
[86]	485	! fin des iles
[4]	486
	487	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
	488	!!$if (itestf.gt.0) then
	489	!!$ write(6,*) 'dans flottab avant dragging asymetrie sur H pour time=',time
	490	!!$ stop
	491	!!$else
	492	!!$ write(6,*) 'dans flottab avant dragging pas d asymetrie sur H pour time=',time
	493	!!$
	494	!!$end if
	495
[86]	496	! 5- calcule les noeuds qui sont streams a l'interieur et donne le betamx et betamy
	497	!----------------------------------------------------------------------------------
[169]	498	call dragging
	499
[86]	500	if (itracebug.eq.1) call tracebug(' routine flottab apres call dragging')
[4]	501	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
	502	!!$if (itestf.gt.0) then
	503	!!$ write(6,*) 'dans flottab apres dragging asymetrie sur H pour time=',time
	504	!!$ stop
	505	!!$else
	506	!!$ write(6,*) 'dans flottab aapres dragging pas d asymetrie sur H pour time=',time
	507	!!$
	508	!!$end if
	509
[86]	510	! 6- calcule les vitesses des points qui sont devenus gzm
	511	!$OMP PARALLEL
	512	!$OMP DO
	513	do j=1,ny
	514	do i=2,nx-1
	515	! si le point etait posé (non gz) et devient gzmx
	516	! definir la direction de la vitesse (moyenne des points)
[4]	517
[86]	518	if ((.not.gz1mx(i,j)).and.(.not.fl1mx(i,j)).and.gzmx(i,j).and. &
	519	(i.gt.2).and.(i.lt.nx)) then
	520	uxs1(i,j)=(uxbar(i+1,j)+uxbar(i-1,j))/2.
	521	endif
[4]	522
[86]	523	end do
	524	end do
	525	!$OMP END DO
[4]	526
[86]	527	!$OMP DO
	528	do j=2,ny-1
	529	do i=1,nx
	530	! si le point etait posé (non gz) et devient gzmy
	531	! definir la direction de la vitesse (moyenne des points)
	532	if ((.not.gz1my(i,j)).and.(.not.fl1my(i,j)).and.gzmy(i,j).and. &
	533	(j.gt.2).and.(j.lt.ny)) then
	534	uys1(i,j)=(uybar(i,j+1)+uybar(i,j-1))/2.
	535	endif
	536
	537	end do
[4]	538	end do
[86]	539	!$OMP END DO
[4]	540
	541
[86]	542	! 7-On determine finalement la position des noeuds stream ou shelf
	543	! -------------------------------------------------------------
[4]	544
[86]	545	if (nt.ge.2) then ! pour ne pas faire ce calcul lors du premier passage
	546	!$OMP WORKSHARE
	547	uxbar(:,:)=uxs1(:,:)
	548	uybar(:,:)=uys1(:,:)
	549	!$OMP END WORKSHARE
	550	endif
[4]	551
[86]	552	!$OMP WORKSHARE
	553	flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:))) &
	554	.or.(.not.marine.and.flotmx(:,:))
[161]	555	where (hmx(:,:).eq.0.)
	556	flgzmx(:,:) = .false.
	557	endwhere
[86]	558	flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:))) &
	559	.or.(.not.marine.and.flotmy(:,:))
[161]	560	where (hmy(:,:).eq.0.)
	561	flgzmy(:,:) = .false.
	562	endwhere
[86]	563	!$OMP END WORKSHARE
[4]	564
	565
[86]	566	! 8- Pour la fusion basale sous les ice shelves- region proche de la grounding line
	567	!---------------------------------------------------------------------------------
	568	! fbm est vrai si le point est flottant mais un des voisins est pose
	569	!_________________________________________________________________________
	570	!$OMP DO
	571	do j=2,ny-1
	572	do i=2,nx-1
[4]	573
[86]	574	! if (i.gt.2.AND.i.lt.nx) then
	575	fbm(i,j)=flot(i,j).and. &
	576	((.not.flot(i+1,j)).or.(.not.flot(i,j+1)) &
	577	.or.(.not.flot(i-1,j)).or.(.not.flot(i,j-1)))
	578	! endif
	579	end do
	580	end do
	581	!$OMP END DO
[4]	582
	583
[86]	584	! 9-On determine maintenant la position du front de glace
	585	! -------------------------------------------------------
	586	! C'est ici que l'on determine la position et le type de front
	587	! print*,'on est dans flottab pour definir les fronts'
[4]	588
	589
[86]	590
[4]	591	!!$do i=3,nx-2
	592	!!$ do j=3,ny-2
	593	!!$ if (h(i,j).gt.1.1) ice(i,j)=1
	594	!!$ end do
	595	!!$end do
[99]	596	! print*, 'flolottab debug', H(71,25),flot(71,25),bm(71,25),ice(71,25),time
	597	!print*,'H(90,179) flottab 2',H(90,179),ice(90,179),flot(90,179)
[86]	598	!$OMP WORKSHARE
	599	where (flot(:,:))
[102]	600	!cdc 1m where (H(:,:).gt.max(Hmin,Hmin+BM(:,:)-Bmelt(:,:)))
[142]	601	where (H(:,:).gt.max(BM(:,:)-Bmelt(:,:)+petit_H,0.)*dt)
	602	!cdc where (H(:,:).gt.0.)
[86]	603	ice(:,:)=1
	604	elsewhere
	605	ice(:,:)=0
[195]	606	H(:,:)=0.
	607	S(:,:)=H(:,:)*(1.-ro/row) + sealevel
	608	B(:,:)=S(:,:) - H(:,:)
[86]	609	end where
	610	elsewhere
[102]	611	where (H(:,:).gt.0.)
[86]	612	ice(:,:)=1
	613	elsewhere
	614	ice(:,:)=0
	615	end where
	616	end where
	617	!$OMP END WORKSHARE
	618	!$OMP END PARALLEL
[195]	619	! print,'flottab',time,H(191,81),bm(191,81),bmelt(191,81),(BM(191,81)-Bmelt(191,81)+petit_H)dt,ice(191,81)
[99]	620	!print,'H(90,179) flottab 3',H(90,179),ice(90,179),flot(90,179),Bsoc(90,179)+H(90,179)ro/row -sealevel
[161]	621
[111]	622	!call determin_front ! cette version ne conserve pas la masse !!!
	623	call determin_front_tof ! version simplifiee
[110]	624
[143]	625	! pour sorties initMIP:
[145]	626	debug_3D(:,:,118) = ice(:,:)*(1-mk(:,:))
[143]	627	debug_3D(:,:,119) = ice(:,:)*mk(:,:)
	628
[99]	629
[86]	630	end subroutine flottab
[4]	631	!--------------------------------------------------------------------
	632
	633	!> SUBROUTINE: determin_tache
	634	!! Routine pour la dtermination du numero de tache a effectuer
	635	!>
	636	subroutine determin_tache
	637
[72]	638	!$ USE OMP_LIB
[71]	639
[4]	640	implicit none
	641	integer :: indice
	642	integer :: label ! no des taches rencontrées dans le mask
	643	integer :: label_max ! no temporaire maxi de tache rencontrées
	644	! integer :: mask_nb = 4
	645	integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales
	646	! integer,dimension(mask_nb) :: mask
	647	integer,dimension(mask_nb) :: mask
	648
	649
	650	! 1-initialisation
	651	!-----------------
	652	label_max=1 ! numero de la tache, la premiere tache est notée 1
	653	label=1
	654	do i=1,n_ta_max
	655	compt(i)=i
	656	enddo
	657	! table_in = .false.
[72]	658	!$OMP PARALLEL
	659	!$OMP WORKSHARE
[4]	660	table_out(:,:) = 0
[111]	661	iceberg1D(:) = .true.
[4]	662	icetrim (:) = .true.
	663	nb_pts_tache(:) = 0
[72]	664	!$OMP END WORKSHARE
	665	!$OMP END PARALLEL
[4]	666	! open(unit=100,file="tache.data",status='replace')
	667
	668	! 2-reperage des taches
	669	!----------------------
	670
[72]	671	do j=2,ny-1
	672	do i=2,nx-1
[4]	673
	674	IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------!
[71]	675
[4]	676	if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then !masque de 2 cases adjacentes
	677	! un des voisins est deja en glace
	678	mask(1) = table_out(i-1,j)
	679	mask(2) = table_out(i,j-1)
	680	label = label_max
	681
	682	!on determine la valeur de la tache minimun (>0) presente ds le masque
	683	do indice=1,mask_nb
	684	if (mask(indice).gt.0) label=min(label,mask(indice))
	685	enddo
[71]	686	!cdc label=min(label,minval(mask(:), mask=mask > 0))
[4]	687
	688	!on fixe la valeur de la tache voisine minimun au point etudie (via label)
	689	table_out(i,j)=label
	690	!si ce noeud est posé, alors la tache n'est pas un iceberg et iceberg=.F.
	691	if (.not.FLOT(I,J)) then
[111]	692	iceberg1D(label)=.false.
[4]	693	endif
	694
	695	!si ce noeud est posé, alors la tache n'est pas un ice stream et icestrim=.F.
	696	if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. &
	697	(.not.gzmy(i,j).and..not.gzmy(i,j+1))) then
	698	icetrim(label)=.false.
	699	endif
	700
	701	! si 2 taches differentes sont dans le masque, il faut les identifier dans compt
[161]	702	! on lui affecte le numero de la tache fondamentale
[4]	703
	704	do indice=1,mask_nb
	705	if(mask(indice).gt.label) then
	706	compt(mask(indice))=-label
	707	endif
	708	enddo
[161]	709	! exemple on est sur le point X : 5 X
	710	do indice=1,mask_nb ! 20
	711	if(mask(indice).gt.label) then ! mask(2)=20 > 5
	712	compt(mask(indice))=label ! compt(20)=5
	713	if (.not.iceberg1D(mask(indice))) iceberg1D(label)=.false. ! si la tache n'etais pas un iceberg => iceberg =.false. iceberg(5)=.false.
	714	if (.not.icetrim(mask(indice))) icetrim(label)=.false.
	715	where (table_out(:,:).eq.mask(indice)) ! where table_out(:,:)=mask(2)=20
	716	table_out(:,:)=label ! table_out(:,:)=label=5
	717	endwhere
	718	endif
	719	enddo
[4]	720
	721	else !aucun des voisins est une tache
	722	table_out(i,j)= label_max
	723	compt(label_max)=label_max
	724	!si ce noeud est posé, alors la ache n'est pas un iceberg et iceberg=.F.
	725	if (.not.FLOT(I,J)) then
[111]	726	iceberg1D(label_max)=.false.
[4]	727	endif
	728
	729	!si ce noeud est posé, alors le tache n'est pas un ice stream et icestrim=.F.
	730	if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. &
	731	(.not.gzmy(i,j).and..not.gzmy(i,j+1))) then
	732	icetrim(label)=.false.
	733	endif
	734
	735	label_max = label_max+1
[161]	736	if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches icebergs=',label_max
[4]	737	endif
	738
	739
	740	else !on est pas sur une tache----------------------------------------------
	741	table_out(i,j)=0 ! Pas necessaire (reecrit 0 sur 0)
	742	endif !---------------------------------------------------------------------
	743
	744
	745	enddo
	746	enddo
	747
	748
[72]	749
[4]	750	! On reorganise compt en ecrivant le numero de la tache fondamentale
	751	! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité)
[71]	752	! On indique aussi le nb de point que contient chaque taches (nb_pts_tache)
	753
[161]	754	!$OMP PARALLEL
	755	!$OMP DO
	756	do j=1,ny
	757	do i=1,nx
[4]	758	if (table_out(i,j).ne.0) then
[161]	759	nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1
[4]	760	endif
[161]	761	enddo
	762	enddo
	763	!$OMP END DO
	764	!$OMP END PARALLEL
[4]	765
[161]	766	!On compte comme englacé uniquement les calottes dont une partie est posée
	767	!$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij)
	768	!$OMP DO
	769	do j=3,ny-2
	770	do i=3,nx-2
	771	test1: if (.not.iceberg1D(table_out(i,j))) then ! on est pas sur un iceberg
	772	if (nb_pts_tache(table_out(i,j)).ge.1) then
	773	ice(i,j)=1
	774	! ici on est sur une tache non iceberg >= 5 points
	775	! on teste si la tache n'est pas completement ice stream
	776
	777	test2: if (icetrim(table_out(i,j))) then ! on a une ile d'ice stream
	778
	779	mask_tache_ij(:,:)=.false.
	780	mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j)) ! pour toute la tache
	781
	782	smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:))
	783	smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:))
	784	smax_i=smax_coord(1)
	785	smax_j=smax_coord(2)
	786
	787	gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false.
	788	gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false.
	789	flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false.
	790	flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false.
	791
	792	if (Smax_.le.sealevel) then
	793	write(num_tracebug,*)'Attention, une ile avec la surface sous l eau'
	794	write(num_tracebug,*)'time=',time,' coord:',smax_i,smax_j
	795	end if
	796	endif test2
	797	end if ! endif deplace
	798	!cdc transfere dans calving :
	799	else ! on est sur un iceberg ! test1
	800	iceberg(i,j)=iceberg1D(table_out(i,j))
	801	!~ ice(i,j)=0
	802	!~ h(i,j)=0. !1. afq, we should put everything in calving!
	803	!~ surnet=H(i,j)*(1.-ro/row)
	804	!~ S(i,j)=surnet+sealevel
	805	!~ B(i,j)=S(i,j)-H(i,j)
	806	endif test1
	807	end do
	808	end do
	809	!$OMP END DO
	810	!$OMP END PARALLEL
[4]	811
[155]	812	debug_3D(:,:,121)=real(table_out(:,:))
	813
[4]	814	end subroutine determin_tache
	815	!----------------------------------------------------------------------
	816	!> SUBROUTINE: determin_front
	817	!!Routine pour la determination du front
	818	!>
	819	subroutine determin_front
[71]	820	!!$ USE OMP_LIB
[4]	821	integer :: i_moins1,i_plus1,i_plus2
	822	integer :: j_moins1,j_plus1,j_plus2
	823
[72]	824	!$OMP PARALLEL
	825	!$OMP DO
[71]	826	do j=3,ny-2
	827	do i=3,nx-2
[4]	828
	829	surice:if (ice(i,j).eq.0) then
	830	do pmx=-1,1,2
	831	do pmy=-1,1,2
	832
	833	diagice : if (ice(i+pmx,j+pmy).eq.1) then
	834
	835	if ((ice(i+pmx,j)+ice(i,j+pmy).eq.2)) then ! test (i) pour eviter les langues
	836	! de glaces diagonales en coin(26dec04)
	837	if ((ice(i+2pmx,j).eq.1.and.ice(i+2pmx,j+pmy).eq.0).or.&
	838	(ice(i,j+2pmy).eq.1.and.ice(i+pmx,j+2pmy).eq.0)) then
	839	ice(i,j)=1
	840	h(i,j)=max(1.,h(i,j))
	841	endif
	842
	843	! test (i) pour eviter les langues de glaces diagonales :
	844	! mouvement du cheval aux echecs
	845
	846	if ((ice(i+2pmx,j+pmy)+ice(i+pmx,j+2pmy).eq.1)) then
	847	if (ice(i+2*pmx,j+pmy).eq.1.and. &
	848	(ice(i+2pmx,j+2pmy)+ice(i,j+2*pmy)).ge.1) then
	849	ice(i,j)=1
	850	h(i,j)=max(1.,h(i,j))
	851	endif
	852	if (ice(i+pmx,j+2*pmy).eq.1.and. &
	853	(ice(i+2pmx,j+2pmy)+ice(i+2*pmx,j)).ge.1) then
	854	ice(i,j)=1
	855	h(i,j)=max(1.,h(i,j))
	856	endif
	857
	858	! test (ii) pour eviter les langues de glaces diagonales :
	859	! le point glace ice(i+pmx,j+pmy) a :
	860	! - ses 4 voisins frontaux en glace
	861	! - mais 2 voisins vides diagonalement opposes
	862
	863	elseif ((ice(i+2pmx,j+pmy)+ice(i+pmx,j+2pmy).eq.2) &
	864	.and.ice(i+2pmx,j+2pmy).eq.0) then
	865
	866	! test (iii) pour faire les tests (i) et (ii)
	867	ice(i,j)=1
	868	h(i,j)=max(1.,h(i,j))
	869	ice(i+2pmx,j+2pmy)=1
	870	h(i+2pmx,j+2pmy)=max(1.,h(i+2pmx,j+2pmy))
	871	endif
	872	endif
	873	endif diagice
	874	enddo
	875	enddo
	876	endif surice
	877	end do
	878	end do
[72]	879	!$OMP END DO
	880	!$OMP ENd PARALLEL
[4]	881
	882	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
	883	!!$if (itestf.gt.0) then
	884	!!$ write(6,*) 'dans front avant remplissage baies asymetrie sur H pour time=',time
	885	!!$ stop
	886	!!$else
	887	!!$ write(6,*) 'dans front avant remplissage baies pas d asymetrie sur H pour time=',time
	888	!!$
	889	!!$end if
	890
	891
	892	! print*,'dans remplissage baies',time
[71]	893
[4]	894	baies: do k=1,2
[72]	895	!$OMP PARALLEL
	896	!$OMP DO PRIVATE(i_moins1,j_moins1,i_plus1,j_plus1,i_plus2,j_plus2)
[4]	897	do j=1,ny
	898	do i=1,nx
	899
	900	surice_xy: if (ice(i,j).eq.0) then
	901	i_moins1=max(i-1,1)
	902	j_moins1=max(j-1,1)
	903	i_plus1=min(i+1,nx)
	904	j_plus1=min(j+1,ny)
	905	i_plus2=min(i+2,nx)
	906	j_plus2=min(j+2,ny)
	907
	908	! test (iii) pour trouver les baies de largeur 1 ou 2 cases
	909	! et combler les trous si ce sont des baies
	910	! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles
	911	! baies horizontales
	912
	913	if (ice(i_moins1,j).eq.1.and.(ice(i_plus1,j).eq.1.or.ice(i_plus2,j).eq.1)) then
	914	if (ice(i,j_moins1).eq.1.or.ice(i,j_plus1).eq.1) then ! ice(i,j)=1
	915	ice(i,j)=1
	916	H(i,j)=max(1.,H(i,j))
	917	endif
	918	endif
	919
	920
	921	if (ice(i,j_moins1).eq.1.and.(ice(i,j_plus1).eq.1.or.ice(i,j_plus2).eq.1)) then
	922	if (ice(i_moins1,j).eq.1.or.ice(i_plus1,j).eq.1) then !ice(i,j)=1
	923	ice(i,j)=1
	924	H(i,j)=max(1.,H(i,j))
	925	endif
	926	endif
	927
	928	endif surice_xy
	929	end do
	930	end do
[72]	931	!$OMP END DO
	932	!$OMP END PARALLEL
[4]	933	end do baies
	934
	935	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
	936	!!$if (itestf.gt.0) then
	937	!!$ write(6,*) 'dans front apres remplissage baies asymetrie sur H pour time=',time
	938	!!$ stop
	939	!!$else
	940	!!$ write(6,*) 'dans front apres remplissage baies pas d asymetrie sur H pour time=',time
	941	!!$
	942	!!$end if
	943
[72]	944	!$OMP PARALLEL
	945	!$OMP DO
[71]	946	do j=2,ny-1
	947	do i=2,nx-1
[4]	948
	949	if (ice(i,j).eq.1) then ! test si ice=1
	950
	951	! if ice, on determine front...
	952	! ainsi, front=0 sur les zones = 0
	953
	954	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
	955	!front= le nb de faces en contact avec un voisin englacé
	956	endif
	957	end do
	958	end do
[72]	959	!$OMP END DO
[4]	960
	961	! traitement des bords. On considere que l'exterieur n'a pas de glace
	962	! attention ce n'est vrai que sur la grande grille
	963
[72]	964	!$OMP DO PRIVATE(i)
[4]	965	do j=2,ny-1
	966	i=1
	967	front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
	968	i=nx
	969	front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
	970	end do
[72]	971	!$OMP END DO
[4]	972
[72]	973	!$OMP DO PRIVATE(j)
[4]	974	do i=2,nx-1
	975	j=1
	976	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))
	977	j=ny
	978	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
	979	end do
[72]	980	!$OMP END DO
	981
[4]	982	! traitement des coins
	983
	984	front(1,1)=ice(2,1)+ice(2,1)
	985	front(1,ny)=ice(2,ny)+ice(1,ny-1)
	986	front(nx,1)=ice(nx,2)+ice(nx-1,1)
	987	front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)
	988
	989	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
	990	!!$if (itestf.gt.0) then
	991	!!$ write(6,*) 'dans front apres front asymetrie sur H pour time=',time
	992	!!$ stop
	993	!!$else
	994	!!$ write(6,*) 'dans front apres front pas d asymetrie sur H pour time=',time
	995	!!$
	996	!!$end if
	997
	998	! on ne compte pas les taches de glace de 2 cases (horizontales ou verticales)
	999	! en fait, si ces deux cases sont flottantes, il faut enlever les icebergs
	1000	! de n'importe quelle taille).
	1001	! si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites
	1002
[72]	1003	!$OMP DO
[4]	1004	do j=1,ny
	1005	do i=1,nx-1
	1006	if (front(i,j).eq.1) then
	1007	if (front(i+1,j).eq.1) then
	1008	ice(i,j)=0
	1009	ice(i+1,j)=0
	1010	front(i,j)=0
	1011	front(i+1,j)=0
	1012	endif
	1013	endif
	1014	end do
	1015	end do
[72]	1016	!$OMP END DO
[4]	1017
[72]	1018	!$OMP DO
[4]	1019	do j=1,ny-1
	1020	do i=1,nx
	1021	if (front(i,j).eq.1) then
	1022	if (front(i,j+1).eq.1) then
	1023	ice(i,j)=0
	1024	ice(i,j+1)=0
	1025	front(i,j)=0
	1026	front(i,j+1)=0
	1027	endif
	1028	end if
	1029	end do
	1030	end do
[72]	1031	!$OMP END DO
[4]	1032
	1033	!isolx signifie pas de voisins en x
	1034	!isoly signifie pas de voisins en y
	1035	!remarque :
	1036	!si isolx/y=.true. alors frontfacex/y=0 (a la fois +1 & -1 or +1-1=0)
	1037
	1038	! calcul de frontfacex et isolx
[72]	1039	!$OMP DO
[4]	1040	do j=1,ny
	1041	do i=2,nx-1
	1042
	1043	if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3)
	1044
	1045	if ((ice(i-1,j)+ice(i+1,j)).lt.2) then ! il y a un front // a x
	1046
	1047	if ((ice(i-1,j)+ice(i+1,j)).eq.0) then
	1048	isolx(i,j)=.true.
	1049	elseif (ice(i-1,j).eq.0) then
	1050	frontfacex(i,j)=-1 ! front i-1 \|i i+1
	1051	else
	1052	frontfacex(i,j)=+1 ! front i-1 i\| i+1
	1053	endif
	1054	endif
	1055	end if !fin du test il y a un front
	1056
	1057	end do
	1058	end do
[72]	1059	!$OMP END DO
[4]	1060
	1061	! calcul de frontfacey et isoly
[72]	1062	!$OMP DO
[4]	1063	do j=2,ny-1
	1064	do i=1,nx
	1065
	1066	if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3)
	1067
	1068	if ((ice(i,j-1)+ice(i,j+1)).lt.2) then ! il y a un front // a y
	1069
	1070	if ((ice(i,j-1)+ice(i,j+1)).eq.0) then
	1071	isoly(i,j)=.true. !front j-1 \|j\| j+1
	1072	elseif (ice(i,j-1).eq.0) then
	1073	frontfacey(i,j)=-1 !front j-1 \|j j+1
	1074	else
	1075	frontfacey(i,j)=+1 !front j-1 j\| j+1
	1076	endif
	1077	endif
	1078	end if !fin du test il y a un front
	1079
	1080	end do
	1081	end do
[72]	1082	!$OMP END DO
[4]	1083
	1084
	1085	! traitement des bords. On considere que l'exterieur n'a pas de glace
	1086	! attention ce n'est vrai que sur la grande grille
	1087
[72]	1088	!$OMP DO PRIVATE(i)
[4]	1089	do j=2,ny-1
	1090	i=1
	1091	if (front(i,j).ge.1) then
	1092	if (ice(i+1,j).eq.0) then
	1093	isolx(i,j)=.true.
	1094	else
	1095	frontfacex(i,j)=-1
	1096	endif
	1097	end if
	1098	i=nx
	1099	if (front(i,j).ge.1) then
	1100	if (ice(i-1,j).eq.0) then
	1101	isolx(i,j)=.true.
	1102	else
	1103	frontfacex(i,j)=1
	1104	endif
	1105	end if
	1106	end do
[72]	1107	!$OMP END DO
[4]	1108
[72]	1109	!$OMP DO PRIVATE(j)
[4]	1110	do i=2,nx-1
	1111	j=1
	1112	if (front(i,j).ge.1) then
	1113	if (ice(i,j+1).eq.0) then
	1114	isoly(i,j)=.true.
	1115	else
	1116	frontfacey(i,j)=-1
	1117	endif
	1118	end if
	1119	j=ny
	1120	if (front(i,j).ge.1) then
	1121	if (ice(i,j-1).eq.0) then
	1122	isoly(i,j)=.true.
	1123	else
	1124	frontfacey(i,j)=1
	1125	endif
	1126	end if
	1127	end do
[72]	1128	!$OMP END DO
	1129	!$OMP END PARALLEL
[4]	1130
	1131	return
	1132	end subroutine determin_front
	1133	!------------------------------------------------------------------------------
	1134
[99]	1135	subroutine determin_front_tof
	1136
	1137	integer,dimension(nx,ny) :: nofront ! tableau de travail (points au dela du front)
	1138	!$OMP PARALLEL
	1139	!$OMP DO
	1140	do j=2,ny-1
	1141	do i=2,nx-1
	1142
	1143	if (ice(i,j).eq.1) then ! test si ice=1
	1144
	1145	! if ice, on determine front...
	1146	! ainsi, front=0 sur les zones = 0
	1147
	1148	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
	1149	!front= le nb de faces en contact avec un voisin englacé
	1150	endif
	1151	end do
	1152	end do
	1153	!$OMP END DO
	1154
	1155	! traitement des bords. On considere que l'exterieur n'a pas de glace
	1156	! attention ce n'est vrai que sur la grande grille
	1157
	1158	!$OMP DO PRIVATE(i)
	1159	do j=2,ny-1
	1160	i=1
	1161	front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
	1162	i=nx
	1163	front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
	1164	end do
	1165	!$OMP END DO
	1166
	1167	!$OMP DO PRIVATE(j)
	1168	do i=2,nx-1
	1169	j=1
	1170	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))
	1171	j=ny
	1172	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
	1173	end do
	1174	!$OMP END DO
[129]	1175	!$OMP BARRIER
[99]	1176	! traitement des coins
	1177
	1178	front(1,1)=ice(2,1)+ice(2,1)
	1179	front(1,ny)=ice(2,ny)+ice(1,ny-1)
	1180	front(nx,1)=ice(nx,2)+ice(nx-1,1)
	1181	front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)
	1182
	1183	! Les points à plus d'un point de grille du bord sont front=-1
[129]	1184	!$OMP WORKSHARE
[99]	1185	nofront(:,:)=0
[129]	1186	!$OMP END WORKSHARE
	1187	!$OMP BARRIER
[99]	1188
	1189	!$OMP DO
	1190	do j=2,ny-1
	1191	do i=2,nx-1
	1192	if ((ice(i,j).eq.0).and.(front(i-1,j)+front(i+1,j)+front(i,j-1)+front(i,j+1)).eq.0) then
	1193	nofront(i,j)=-1
	1194	endif
	1195	enddo
	1196	enddo
	1197	!$OMP END DO
[129]	1198	!$OMP BARRIER
	1199
	1200	!$OMP WORKSHARE
[99]	1201	where (nofront(:,:).eq.-1)
	1202	front(:,:)=-1
	1203	endwhere
[129]	1204	!$OMP END WORKSHARE
[99]	1205	!$OMP END PARALLEL
	1206
	1207	end subroutine determin_front_tof
	1208
[110]	1209
	1210	!> SUBROUTINE: determin_marais
	1211	!! afq -- Routine pour l'identification des "marais"
	1212	!! Un marais est une tache "shelf" entouré de points grounded
	1213	!! Copie sauvage de determin_tache, adapte au probleme du marais
	1214	!>
	1215	subroutine determin_marais
	1216
[156]	1217	!$ USE OMP_LIB
[110]	1218
[156]	1219	implicit none
	1220
	1221	integer :: indice
	1222	integer :: label ! no des taches rencontrées dans le mask
	1223	integer :: label_max ! no temporaire maxi de tache rencontrées
	1224	integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales
	1225	integer,dimension(mask_nb) :: mask ! numero de tache des points adjacents
[110]	1226
[156]	1227	integer,dimension(nx,ny) :: table_out_marais !< numeros de tache d'un point ij
	1228	integer,dimension(0:n_ta_max) :: compt_marais !< contient les equivalence entre les taches
	1229	integer,dimension(0:n_ta_max) :: nb_pts_marais !< indique le nombre de points par tache
	1230	logical,dimension(0:n_ta_max) :: marais !< T si flottants entoure de poses, F sinon
	1231
	1232
[110]	1233	! 1-initialisation
	1234	!-----------------
[156]	1235	label_max=1 ! numero de la tache, la premiere tache est notée 1
	1236	label=1
	1237	do i=1,n_ta_max
	1238	compt_marais(i)=i
	1239	enddo
	1240	!$OMP PARALLEL
	1241	!$OMP WORKSHARE
	1242	table_out_marais(:,:) = 0
	1243	marais(:) = .true.
	1244	nb_pts_marais(:) = 0
	1245	!$OMP END WORKSHARE
	1246	!$OMP END PARALLEL
[110]	1247
	1248	! 2-reperage des taches
	1249	!----------------------
	1250
[156]	1251	do j=2,ny-1
	1252	do i=2,nx-1
	1253	if ((ice(i,j).ge.1).and.flot(i,j)) then ! on est sur la glace qui flotte-------------------!
[110]	1254
	1255	if (((ice(i-1,j).ge.1).and.flot(i-1,j)).or.((ice(i,j-1).ge.1).and.flot(i,j-1))) then !masque de 2 cases adjacentes
[156]	1256	! un des voisins est deja en glace
	1257	mask(1) = table_out_marais(i-1,j)
	1258	mask(2) = table_out_marais(i,j-1)
	1259	label = label_max
[110]	1260
[156]	1261	! on determine la valeur de la tache minimun (>0) presente ds le masque
	1262	do indice=1,mask_nb
	1263	if (mask(indice).gt.0) label=min(label,mask(indice))
	1264	enddo
[110]	1265
[156]	1266	! on fixe la valeur de la tache voisine minimun au point etudie (via label)
	1267	table_out_marais(i,j)=label
[110]	1268
[156]	1269	!si un des voisins n'est pas glace alors la tache n'est pas un marais
	1270	if ( (ice(i+1,j).eq.0) .or. (ice(i,j+1).eq.0) .or. (ice(i-1,j).eq.0) .or. (ice(i,j-1).eq.0) ) then
	1271	marais(label)=.false.
	1272	endif
[110]	1273
[156]	1274	! si 2 taches differentes sont dans le masque, il faut les identifier dans compt_marais
	1275	! on lui affecte le numero de la tache fondamentale
[110]	1276
[156]	1277	! exemple on est sur le point X : 5 X
	1278	do indice=1,mask_nb ! 20
	1279	if(mask(indice).gt.label) then ! mask(2)=20 > 5
	1280	compt_marais(mask(indice))=label ! compt_marais(20)=5
	1281	if (.not.marais(mask(indice))) marais(label)=.false. ! si la tache n'etais pas un marais => marais =.false. marais(-(-5))=.false.
	1282	where (table_out_marais(:,:).eq.mask(indice)) ! where table_out_marais(:,:)=mask(2)=20
	1283	table_out_marais(:,:)=label ! table_out_marais(:,:)=label=5
	1284	endwhere
	1285	endif
	1286	enddo
[110]	1287
	1288	else !aucun des voisins est une tache
[156]	1289	table_out_marais(i,j)= label_max
	1290	compt_marais(label_max)=label_max
[110]	1291
[156]	1292	! si un des voisins n'est pas glace alors la tache n'est pas un marais
	1293	if ( (ice(i+1,j).eq.0) .or. (ice(i,j+1).eq.0) .or. (ice(i-1,j).eq.0) .or. (ice(i,j-1).eq.0) ) then
	1294	marais(label_max)=.false.
	1295	endif
	1296	label_max = label_max+1
	1297	if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches=',label_max
[110]	1298	endif
[156]	1299	else ! on est pas sur une tache--------------------------------------------
[110]	1300	table_out_marais(i,j)=0 ! Pas necessaire (reecrit 0 sur 0)
[156]	1301	endif !---------------------------------------------------------------------
	1302	enddo
[110]	1303	enddo
	1304
[156]	1305	marais(0)=.false.
[110]	1306
[156]	1307	!$OMP PARALLEL
	1308	!$OMP DO
	1309	do j=1,ny
	1310	do i=1,nx
[110]	1311	if (table_out_marais(i,j).ne.0) then
[156]	1312	nb_pts_marais(compt_marais(table_out_marais(i,j)))= nb_pts_marais(compt_marais(table_out_marais(i,j)))+1
[110]	1313	endif
	1314	flot_marais(i,j) = marais(table_out_marais(i,j))
[156]	1315	enddo
	1316	enddo
	1317	!$OMP END DO
	1318	!$OMP END PARALLEL
[110]	1319
[156]	1320	debug_3D(:,:,122)=real(table_out_marais(:,:))
[154]	1321
[110]	1322	end subroutine determin_marais
	1323
[4]	1324	end module flottab_mod
	1325
	1326

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

Download in other formats: