Context Navigation

source: branches/GRISLIv3/SOURCES/flottab2-0.7.f90 @ 406

Last change on this file since 406 was 396, checked in by dumas, 16 months ago
use only in module flottab_mod
File size: 46.3 KB

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

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

Download in other formats: