Context Navigation

flottab2-0.7.f90 @ 334

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

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

Context Navigation

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

Download in other formats: