Context Navigation

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

Last change on this file since 467 was 467, checked in by aquiquet, 4 months ago
Cleaning branch: continuing module3D cleaning
File size: 43.6 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 geography, 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,&
107	ice,front,&
108	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	logical,dimension(nx,ny) :: new_flot_point !< pour signaler les points qui se mettent a flotter entre 2 pas de temps dtt
116	logical,dimension(nx,ny) :: new_flotmx !< pour signaler les points qui deviennent flottantmx entre 2 dtt
117	logical,dimension(nx,ny) :: new_flotmy !< pour signaler les points qui deviennent flottantmy entre 2 dtt
118
119	integer :: i,j
120
121	if (itracebug.eq.1) call tracebug(' Entree dans routine flottab')
122
123	SHELFY = .FALSE.
124
125
126	! cas particulier des runs paleo ou on impose un masque grounded
127
128	!$OMP PARALLEL PRIVATE(archim,surnet)
129	if (igrdline.eq.2) then
130	!$OMP WORKSHARE
131	where ( mk_init(:,:).eq.1) ! pose
132	flot(:,:) = .False.
133	H(:,:)=max(H(:,:),(10.+sealevel_2d(:,:)-Bsoc(:,:))*row/ro) ! pour avoir archim=10
134	S(:,:) = Bsoc(:,:) + H(:,:)
135	elsewhere
136	flot(:,:) = .True.
137	end where
138	!$OMP END WORKSHARE
139	end if
140
141
142	! 1-INITIALISATION
143	! ----------------
144	! initialisation des variables pour detecter les points qui se mettent
145	! a flotter entre 2 dtt
146
147	!$OMP DO
148	do j=1,ny
149	do i=1,nx
150	new_flot_point(i,j)=.false.
151	new_flotmx(i,j)=.false.
152	new_flotmy(i,j)=.false.
153	enddo
154	enddo
155	!$OMP END DO
156
157	! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m
158
159	!$OMP WORKSHARE
160	ICE(:,:)=0
161	front(:,:)=0
162	iceberg(:,:)=.false.
163	!$OMP END WORKSHARE
164
165	! fin de l'initialisation
166	!_____________________________________________________________________
167
168	! 2-TESTE LES NOUVEAUX POINTS FLOTTANTS
169	! -------------------------------------
170
171	!$OMP DO
172	do j=1,ny
173	do i=1,nx
174
175	uxs1(i,j)=uxbar(i,j)
176	uys1(i,j)=uybar(i,j)
177
178	archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j)
179	! if ((i.eq.132).and.(j.eq.183)) print*,'archim=',archim
180
181
182	arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.e-3)) then ! le point flotte
183	mk(i,j)=1
184
185
186	ex_pose: if ((.not.FLOT(I,J)).and.(isynchro.eq.1)) then ! il ne flottait pas avant
187	FLOT(I,J)=.true.
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	endif
218	!cdc correction topo pour suivre variations sealevel
219	!cdd S(i,j)=Bsoc(i,j)+H(i,j)
220	S(i,j)=Bsoc(i,j)+H(i,j)
221	B(i,j)=Bsoc(i,j)
222
223	else if ((H(i,j).LE.0.).and.(archim.LT.0.)) then ! terre deglace qui devient ocean
224	!cdc ice(i,j)=0
225	!cdc H(i,j)=1.
226	!cdc 1m H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01))
227	flot(i,j)=.true. !cdc points ocean sont flot meme sans glace
228	H(i,j)=0.
229	S(i,j)=sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
230	B(i,j)=S(i,j)-H(i,j)
231	endif arch
232
233	! Si la glace flotte -> PRUDENCE !!!
234	! S et B sont alors determines avec la condition de flottabilite
235
236	if (flot(i,j)) then
237	shelfy = .true.
238
239	surnet=H(i,j)*(1.-ro/row)
240	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?
241	B(i,j)=S(i,j)-H(i,j)
242	end if
243
244	end do
245	end do
246	!$OMP END DO
247
248	!!$ do i=1,nx
249	!!$ do j=1,ny
250	!!$ if (flot(i,j)) then
251	!!$ mk(i,j)=1
252	!!$ else
253	!!$ mk(i,j)=0
254	!!$ endif
255	!!$ end do
256	!!$ end do
257
258
259
260	!-----------------------------------------------------------------------
261	!$OMP DO
262	domain_x: do j=1,ny
263	do i=2,nx
264
265	! 3_x A- NOUVELLE DEFINITION DE FLOTMX, LES POINTS
266	! AYANT UN DES VOISINS FLOTTANTS SONT FLOTMX ET GZMX
267	! -------------------------------------------
268
269	! fl1 est l'ancienne valeur de flotmx
270	gz1mx(i,j)=gzmx(i,j) ! gz1 est l'ancienne valeur de gzmx
271	fl1mx(i,j)=flotmx(i,j) ! fl1 est l'ancienne valeur de flotmx
272
273	flotmx(i,j)=flot(i,j).and.flot(i-1,j)
274
275	! test pour detecter les nouveaux flotmx entre 2 dtt :
276
277	if (flotmx(i,j).and.(new_flot_point(i,j).or. &
278	new_flot_point(i-1,j))) then
279	new_flotmx(i,j)=.true.
280	endif
281
282
283	! premiere determination de gzmx
284	!__________________________________________________________________________
285
286	! gzmx si un des deux voisins est flottant et l'autre posé
287	! i-1 i
288	gzmx(i,j)=((flot(i,j).and..not.flot(i-1,j)) & ! F P
289	.or.(.not.flot(i,j).and.flot(i-1,j))) ! P F
290
291	! A condition d'etre assez proche de la flottaison
292	! sur le demi noeud condition archim < 100 m
293
294	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)
295	gzmx(i,j)=gzmx(i,j).and.(archim.le.100.)
296
297	end do
298	end do domain_x
299	!$OMP END DO
300	!if (itracebug.eq.1) call tracebug(' routine flottab apres domain_x')
301
302	! 3_y B- NOUVELLE DEFINITION DE FLOTMY
303	! --------------------------------
304	!$OMP DO
305	domain_y: do j=2,ny
306	do i=1,nx
307
308	gz1my(i,j)=gzmy(i,j) ! gz1 est l'ancienne valeur de gzmy
309	fl1my(i,j)=flotmy(i,j) ! fl1 est l'ancienne valeur de flotmy
310
311	flotmy(i,j)=flot(i,j).and.flot(i,j-1)
312
313	! test pour detecter les nouveaux flotmy entre 2 dtt :
314
315	if (flotmy(i,j).and.(new_flot_point(i,j).or. &
316	new_flot_point(i,j-1))) then
317	new_flotmy(i,j)=.true.
318	endif
319
320	! premiere determination de gzmy
321	!__________________________________________________________________________
322
323	! gzmy si un des deux voisins est flottant et l'autre posé
324
325	gzmy(i,j)=((flot(i,j).and..not.flot(i,j-1)) &
326	.or.(.not.flot(i,j).and.flot(i,j-1)))
327
328	! A condition d'etre assez proche de la flottaison
329	! sur le demi noeud condition archim > 100 m
330
331	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)
332	gzmy(i,j)=gzmy(i,j).and.(archim.le.100.)
333
334	end do
335	end do domain_y
336	!$OMP END DO
337
338
339	!-------------------------------------------------------------------------------------
340	! attention : pour expériences Heino
341	! gzmy(i,j)=gzmy_heino(i,j)
342	! shelfy=.true.
343	! shelfy=.false.
344	!____________________________________________________________________________________
345
346
347	!!$
348	!!$
349	!!$! 4- Condition sur les bords
350	!!$
351	!!$
352	!!$ do i=2,nx
353	!!$ flotmx(i,1) = (flot(i,1).or.flot(i-1,1))
354	!!$ flotmy(i,1) = .false.
355	!!$ end do
356	!!$
357	!!$ do j=2,ny
358	!!$ flotmy(1,j) = (flot(1,j).or.flot(1,j-1))
359	!!$ flotmx(1,j) = .false.
360	!!$ end do
361	!!$
362	!!$ flotmx(1,1) = .false.
363	!!$ flotmy(1,1) = .false.
364	!!$
365
366
367	! 4- determination des iles
368	! -------------------------
369	!$OMP WORKSHARE
370	ilemx(:,:)=.false.
371	ilemy(:,:)=.false.
372	!$OMP END WORKSHARE
373
374	! afq -- 17/01/19: on supprime les iles.
375	! ! selon x
376	! !$OMP DO
377	! ilesx: do j=2,ny-1
378	! do i=3,nx-2
379	! ! F G F
380	! ! x
381	! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
382	! if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. &
383	! (sdx(i,j).LT.1.E-02)) then
384	! ilemx(i,j)=.true.
385	! ilemx(i+1,j)=.true.
386
387	! ! F G G F
388	! ! x
389	! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
390	! else if ((flot(i-1,j).and..not.flot(i,j) &
391	! .and..not.flot(i+1,j)).and.flot(i+2,j).and. &
392	! (sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then
393	! ilemx(i,j)=.true.
394	! ilemx(i+1,j)=.true.
395	! ilemx(i+2,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-2,j).and..not.flot(i-1,j) &
401	! .and..not.flot(i,j)).and.flot(i+1,j).and. &
402	! (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then
403	! ilemx(i-1,j)=.true.
404	! ilemx(i,j)=.true.
405	! ilemx(i+1,j)=.true.
406
407	! ! F G G G F
408	! ! x
409	! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
410	! else if ((i.lt.nx-2) &
411	! .and.(flot(i-2,j).and..not.flot(i-1,j) &
412	! .and..not.flot(i,j)).and..not.flot(i+1,j) &
413	! .and.flot(i+2,j).and. &
414	! (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02 &
415	! .and.sdx(i+1,j).LT.1.E-02)) then
416	! ilemx(i-1,j)=.true.
417	! ilemx(i,j)=.true.
418	! ilemx(i+1,j)=.true.
419	! ilemx(i+2,j)=.true.
420
421	! endif
422
423	! end do
424	! end do ilesx
425	! !$OMP END DO
426
427	! ! selon y
428	! !$OMP DO
429	! ilesy: do j=3,ny-2
430	! do i=2,nx-1
431	! ! F G F
432	! ! x
433	! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
434	! if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. &
435	! (sdy(i,j).LT.1.E-02)) then
436	! ilemy(i,j)=.true.
437	! ilemy(i,j+1)=.true.
438
439	! ! F G G F
440	! ! x
441	! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
442	! else if ((flot(i,j-1).and..not.flot(i,j) &
443	! .and..not.flot(i,j+1)).and.flot(i,j+2).and. &
444	! (sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then
445	! ilemy(i,j)=.true.
446	! ilemy(i,j+1)=.true.
447	! ilemy(i,j+2)=.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-2).and..not.flot(i,j-1) &
453	! .and..not.flot(i,j)).and.flot(i,j+1).and. &
454	! (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then
455	! ilemy(i,j-1)=.true.
456	! ilemy(i,j)=.true.
457	! ilemy(i,j+1)=.true.
458
459	! ! F G G G F
460	! ! x
461	! ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
462	! else if ((j.lt.ny-2) &
463	! .and.(flot(i,j-2).and..not.flot(i,j-1) &
464	! .and..not.flot(i,j)).and..not.flot(i,j+1) &
465	! .and.flot(i,j+2).and. &
466	! (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02 &
467	! .and.sdy(i,j+1).LT.1.E-02)) then
468	! ilemy(i,j-1)=.true.
469	! ilemy(i,j)=.true.
470	! ilemy(i,j+1)=.true.
471	! ilemy(i,j+2)=.true.
472	! endif
473	! end do
474	! end do ilesy
475	! !$OMP END DO
476	! ! fin des iles
477
478	!$OMP END PARALLEL
479
480
481	! 5- calcule les noeuds qui sont streams a l'interieur et donne le betamx et betamy
482	!----------------------------------------------------------------------------------
483	call mstream_dragging
484
485	if (itracebug.eq.1) call tracebug(' routine flottab apres call dragging')
486	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
487	!!$if (itestf.gt.0) then
488	!!$ write(6,*) 'dans flottab apres dragging asymetrie sur H pour time=',time
489	!!$ stop
490	!!$else
491	!!$ write(6,*) 'dans flottab aapres dragging pas d asymetrie sur H pour time=',time
492	!!$
493	!!$end if
494
495	! 6- calcule les vitesses des points qui sont devenus gzm
496	!$OMP PARALLEL
497	!$OMP DO
498	do j=1,ny
499	do i=2,nx-1
500	! si le point etait posé (non gz) et devient gzmx
501	! definir la direction de la vitesse (moyenne des points)
502
503	if ((.not.gz1mx(i,j)).and.(.not.fl1mx(i,j)).and.gzmx(i,j).and. &
504	(i.gt.2).and.(i.lt.nx)) then
505	uxs1(i,j)=(uxbar(i+1,j)+uxbar(i-1,j))/2.
506	endif
507
508	end do
509	end do
510	!$OMP END DO
511
512	!$OMP DO
513	do j=2,ny-1
514	do i=1,nx
515	! si le point etait posé (non gz) et devient gzmy
516	! definir la direction de la vitesse (moyenne des points)
517	if ((.not.gz1my(i,j)).and.(.not.fl1my(i,j)).and.gzmy(i,j).and. &
518	(j.gt.2).and.(j.lt.ny)) then
519	uys1(i,j)=(uybar(i,j+1)+uybar(i,j-1))/2.
520	endif
521
522	end do
523	end do
524	!$OMP END DO
525
526
527	! 7-On determine finalement la position des noeuds stream ou shelf
528	! -------------------------------------------------------------
529
530	if (nt.ge.2) then ! pour ne pas faire ce calcul lors du premier passage
531	!$OMP WORKSHARE
532	uxbar(:,:)=uxs1(:,:)
533	uybar(:,:)=uys1(:,:)
534	!$OMP END WORKSHARE
535	endif
536
537	!$OMP WORKSHARE
538	flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:))) &
539	.or.(.not.marine.and.flotmx(:,:))
540	where (hmx(:,:).eq.0.)
541	flgzmx(:,:) = .false.
542	endwhere
543	flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:))) &
544	.or.(.not.marine.and.flotmy(:,:))
545	where (hmy(:,:).eq.0.)
546	flgzmy(:,:) = .false.
547	endwhere
548	!$OMP END WORKSHARE
549
550
551	! 8- Pour la fusion basale sous les ice shelves- region proche de la grounding line
552	!---------------------------------------------------------------------------------
553	! fbm est vrai si le point est flottant mais un des voisins est pose
554	!_________________________________________________________________________
555	!$OMP DO
556	do j=2,ny-1
557	do i=2,nx-1
558
559	! if (i.gt.2.AND.i.lt.nx) then
560	fbm(i,j)=flot(i,j).and. &
561	((.not.flot(i+1,j)).or.(.not.flot(i,j+1)) &
562	.or.(.not.flot(i-1,j)).or.(.not.flot(i,j-1)))
563	! endif
564	end do
565	end do
566	!$OMP END DO
567
568
569	! 9-On determine maintenant la position du front de glace
570	! -------------------------------------------------------
571	! C'est ici que l'on determine la position et le type de front
572	! print*,'on est dans flottab pour definir les fronts'
573
574
575
576	!!$do i=3,nx-2
577	!!$ do j=3,ny-2
578	!!$ if (h(i,j).gt.1.1) ice(i,j)=1
579	!!$ end do
580	!!$end do
581	! print*, 'flolottab debug', H(71,25),flot(71,25),bm(71,25),ice(71,25),time
582	!print*,'H(90,179) flottab 2',H(90,179),ice(90,179),flot(90,179)
583	!$OMP WORKSHARE
584	where (flot(:,:))
585	!cdc 1m where (H(:,:).gt.max(Hmin,Hmin+BM(:,:)-Bmelt(:,:)))
586	where (H(:,:).gt.max(BM(:,:)-Bmelt(:,:)+petit_H,0.)*dt)
587	!cdc where (H(:,:).gt.0.)
588	ice(:,:)=1
589	elsewhere
590	ice(:,:)=0
591	H(:,:)=0.
592	S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:)
593	B(:,:)=S(:,:) - H(:,:)
594	end where
595	elsewhere
596	where (H(:,:).gt.0.)
597	ice(:,:)=1
598	elsewhere
599	ice(:,:)=0
600	end where
601	end where
602	!$OMP END WORKSHARE
603	!$OMP END PARALLEL
604	! print,'flottab',time,H(191,81),bm(191,81),bmelt(191,81),(BM(191,81)-Bmelt(191,81)+petit_H)dt,ice(191,81)
605	!print,'H(90,179) flottab 3',H(90,179),ice(90,179),flot(90,179),Bsoc(90,179)+H(90,179)ro/row -sealevel
606
607	!call determin_front ! cette version ne conserve pas la masse !!!
608	call determin_front_tof ! version simplifiee
609
610	! pour sorties initMIP:
611	debug_3D(:,:,118) = ice(:,:)*(1-mk(:,:))
612	debug_3D(:,:,119) = ice(:,:)*mk(:,:)
613
614
615	end subroutine flottab
616	!--------------------------------------------------------------------
617
618	!> SUBROUTINE: determin_tache
619	!! Routine pour la dtermination du numero de tache a effectuer
620	!>
621	subroutine determin_tache
622
623	use module3D_phy, only:ice,flot,gzmx,gzmy,S,flgzmx,flgzmy,sealevel,time,iceberg,&
624	debug_3D
625	use runparam, only: num_tracebug
626
627	!$ USE OMP_LIB
628
629	implicit none
630	integer :: i,j
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
812	use module3D_phy, only:ice,H,front
813
814
815	!!$ USE OMP_LIB
816	implicit none
817
818	integer :: i,j,k
819	integer :: i_moins1,i_plus1,i_plus2
820	integer :: j_moins1,j_plus1,j_plus2
821
822	!$OMP PARALLEL
823	!$OMP DO
824	do j=3,ny-2
825	do i=3,nx-2
826
827	surice:if (ice(i,j).eq.0) then
828	do pmx=-1,1,2
829	do pmy=-1,1,2
830
831	diagice : if (ice(i+pmx,j+pmy).eq.1) then
832
833	if ((ice(i+pmx,j)+ice(i,j+pmy).eq.2)) then ! test (i) pour eviter les langues
834	! de glaces diagonales en coin(26dec04)
835	if ((ice(i+2pmx,j).eq.1.and.ice(i+2pmx,j+pmy).eq.0).or.&
836	(ice(i,j+2pmy).eq.1.and.ice(i+pmx,j+2pmy).eq.0)) then
837	ice(i,j)=1
838	h(i,j)=max(1.,h(i,j))
839	endif
840
841	! test (i) pour eviter les langues de glaces diagonales :
842	! mouvement du cheval aux echecs
843
844	if ((ice(i+2pmx,j+pmy)+ice(i+pmx,j+2pmy).eq.1)) then
845	if (ice(i+2*pmx,j+pmy).eq.1.and. &
846	(ice(i+2pmx,j+2pmy)+ice(i,j+2*pmy)).ge.1) then
847	ice(i,j)=1
848	h(i,j)=max(1.,h(i,j))
849	endif
850	if (ice(i+pmx,j+2*pmy).eq.1.and. &
851	(ice(i+2pmx,j+2pmy)+ice(i+2*pmx,j)).ge.1) then
852	ice(i,j)=1
853	h(i,j)=max(1.,h(i,j))
854	endif
855
856	! test (ii) pour eviter les langues de glaces diagonales :
857	! le point glace ice(i+pmx,j+pmy) a :
858	! - ses 4 voisins frontaux en glace
859	! - mais 2 voisins vides diagonalement opposes
860
861	elseif ((ice(i+2pmx,j+pmy)+ice(i+pmx,j+2pmy).eq.2) &
862	.and.ice(i+2pmx,j+2pmy).eq.0) then
863
864	! test (iii) pour faire les tests (i) et (ii)
865	ice(i,j)=1
866	h(i,j)=max(1.,h(i,j))
867	ice(i+2pmx,j+2pmy)=1
868	h(i+2pmx,j+2pmy)=max(1.,h(i+2pmx,j+2pmy))
869	endif
870	endif
871	endif diagice
872	enddo
873	enddo
874	endif surice
875	end do
876	end do
877	!$OMP END DO
878	!$OMP ENd PARALLEL
879
880	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
881	!!$if (itestf.gt.0) then
882	!!$ write(6,*) 'dans front avant remplissage baies asymetrie sur H pour time=',time
883	!!$ stop
884	!!$else
885	!!$ write(6,*) 'dans front avant remplissage baies pas d asymetrie sur H pour time=',time
886	!!$
887	!!$end if
888
889
890	! print*,'dans remplissage baies',time
891
892	baies: do k=1,2
893	!$OMP PARALLEL
894	!$OMP DO PRIVATE(i_moins1,j_moins1,i_plus1,j_plus1,i_plus2,j_plus2)
895	do j=1,ny
896	do i=1,nx
897
898	surice_xy: if (ice(i,j).eq.0) then
899	i_moins1=max(i-1,1)
900	j_moins1=max(j-1,1)
901	i_plus1=min(i+1,nx)
902	j_plus1=min(j+1,ny)
903	i_plus2=min(i+2,nx)
904	j_plus2=min(j+2,ny)
905
906	! test (iii) pour trouver les baies de largeur 1 ou 2 cases
907	! et combler les trous si ce sont des baies
908	! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles
909	! baies horizontales
910
911	if (ice(i_moins1,j).eq.1.and.(ice(i_plus1,j).eq.1.or.ice(i_plus2,j).eq.1)) then
912	if (ice(i,j_moins1).eq.1.or.ice(i,j_plus1).eq.1) then ! ice(i,j)=1
913	ice(i,j)=1
914	H(i,j)=max(1.,H(i,j))
915	endif
916	endif
917
918
919	if (ice(i,j_moins1).eq.1.and.(ice(i,j_plus1).eq.1.or.ice(i,j_plus2).eq.1)) then
920	if (ice(i_moins1,j).eq.1.or.ice(i_plus1,j).eq.1) then !ice(i,j)=1
921	ice(i,j)=1
922	H(i,j)=max(1.,H(i,j))
923	endif
924	endif
925
926	endif surice_xy
927	end do
928	end do
929	!$OMP END DO
930	!$OMP END PARALLEL
931	end do baies
932
933	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
934	!!$if (itestf.gt.0) then
935	!!$ write(6,*) 'dans front apres remplissage baies asymetrie sur H pour time=',time
936	!!$ stop
937	!!$else
938	!!$ write(6,*) 'dans front apres remplissage baies pas d asymetrie sur H pour time=',time
939	!!$
940	!!$end if
941
942	!$OMP PARALLEL
943	!$OMP DO
944	do j=2,ny-1
945	do i=2,nx-1
946
947	if (ice(i,j).eq.1) then ! test si ice=1
948
949	! if ice, on determine front...
950	! ainsi, front=0 sur les zones = 0
951
952	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
953	!front= le nb de faces en contact avec un voisin englacé
954	endif
955	end do
956	end do
957	!$OMP END DO
958
959	! traitement des bords. On considere que l'exterieur n'a pas de glace
960	! attention ce n'est vrai que sur la grande grille
961
962	!$OMP DO PRIVATE(i)
963	do j=2,ny-1
964	i=1
965	front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
966	i=nx
967	front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
968	end do
969	!$OMP END DO
970
971	!$OMP DO PRIVATE(j)
972	do i=2,nx-1
973	j=1
974	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))
975	j=ny
976	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
977	end do
978	!$OMP END DO
979
980	! traitement des coins
981
982	front(1,1)=ice(2,1)+ice(2,1)
983	front(1,ny)=ice(2,ny)+ice(1,ny-1)
984	front(nx,1)=ice(nx,2)+ice(nx-1,1)
985	front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)
986
987	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
988	!!$if (itestf.gt.0) then
989	!!$ write(6,*) 'dans front apres front asymetrie sur H pour time=',time
990	!!$ stop
991	!!$else
992	!!$ write(6,*) 'dans front apres front pas d asymetrie sur H pour time=',time
993	!!$
994	!!$end if
995
996	! on ne compte pas les taches de glace de 2 cases (horizontales ou verticales)
997	! en fait, si ces deux cases sont flottantes, il faut enlever les icebergs
998	! de n'importe quelle taille).
999	! si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites
1000
1001	!$OMP DO
1002	do j=1,ny
1003	do i=1,nx-1
1004	if (front(i,j).eq.1) then
1005	if (front(i+1,j).eq.1) then
1006	ice(i,j)=0
1007	ice(i+1,j)=0
1008	front(i,j)=0
1009	front(i+1,j)=0
1010	endif
1011	endif
1012	end do
1013	end do
1014	!$OMP END DO
1015
1016	!$OMP DO
1017	do j=1,ny-1
1018	do i=1,nx
1019	if (front(i,j).eq.1) then
1020	if (front(i,j+1).eq.1) then
1021	ice(i,j)=0
1022	ice(i,j+1)=0
1023	front(i,j)=0
1024	front(i,j+1)=0
1025	endif
1026	end if
1027	end do
1028	end do
1029	!$OMP END DO
1030
1031	!$OMP END PARALLEL
1032
1033	return
1034	end subroutine determin_front
1035	!------------------------------------------------------------------------------
1036
1037	subroutine determin_front_tof
1038
1039	use module3D_phy, only:front,ice
1040
1041	implicit none
1042
1043	integer :: i,j
1044	integer,dimension(nx,ny) :: nofront ! tableau de travail (points au dela du front)
1045	!$OMP PARALLEL
1046	!$OMP DO
1047	do j=2,ny-1
1048	do i=2,nx-1
1049
1050	if (ice(i,j).eq.1) then ! test si ice=1
1051
1052	! if ice, on determine front...
1053	! ainsi, front=0 sur les zones = 0
1054
1055	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
1056	!front= le nb de faces en contact avec un voisin englacé
1057	endif
1058	end do
1059	end do
1060	!$OMP END DO
1061
1062	! traitement des bords. On considere que l'exterieur n'a pas de glace
1063	! attention ce n'est vrai que sur la grande grille
1064
1065	!$OMP DO PRIVATE(i)
1066	do j=2,ny-1
1067	i=1
1068	front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
1069	i=nx
1070	front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
1071	end do
1072	!$OMP END DO
1073
1074	!$OMP DO PRIVATE(j)
1075	do i=2,nx-1
1076	j=1
1077	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))
1078	j=ny
1079	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
1080	end do
1081	!$OMP END DO
1082	!$OMP BARRIER
1083	! traitement des coins
1084
1085	front(1,1)=ice(2,1)+ice(2,1)
1086	front(1,ny)=ice(2,ny)+ice(1,ny-1)
1087	front(nx,1)=ice(nx,2)+ice(nx-1,1)
1088	front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)
1089
1090	! Les points à plus d'un point de grille du bord sont front=-1
1091	!$OMP WORKSHARE
1092	nofront(:,:)=0
1093	!$OMP END WORKSHARE
1094	!$OMP BARRIER
1095
1096	!$OMP DO
1097	do j=2,ny-1
1098	do i=2,nx-1
1099	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
1100	nofront(i,j)=-1
1101	endif
1102	enddo
1103	enddo
1104	!$OMP END DO
1105	!$OMP BARRIER
1106
1107	!$OMP WORKSHARE
1108	where (nofront(:,:).eq.-1)
1109	front(:,:)=-1
1110	endwhere
1111	!$OMP END WORKSHARE
1112	!$OMP END PARALLEL
1113
1114	end subroutine determin_front_tof
1115
1116
1117	!> SUBROUTINE: determin_marais
1118	!! afq -- Routine pour l'identification des "marais"
1119	!! Un marais est une tache "shelf" entouré de points grounded
1120	!! Copie sauvage de determin_tache, adapte au probleme du marais
1121	!>
1122	subroutine determin_marais
1123
1124	use module3D_phy, only:ice,flot,flot_marais,debug_3D
1125	!$ USE OMP_LIB
1126
1127	implicit none
1128
1129	integer :: i,j
1130	integer :: indice
1131	integer :: label ! no des taches rencontrées dans le mask
1132	integer :: label_max ! no temporaire maxi de tache rencontrées
1133	integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales
1134	integer,dimension(mask_nb) :: mask ! numero de tache des points adjacents
1135
1136	integer,dimension(nx,ny) :: table_out_marais !< numeros de tache d'un point ij
1137	integer,dimension(0:n_ta_max) :: compt_marais !< contient les equivalence entre les taches
1138	integer,dimension(0:n_ta_max) :: nb_pts_marais !< indique le nombre de points par tache
1139	logical,dimension(0:n_ta_max) :: marais !< T si flottants entoure de poses, F sinon
1140
1141
1142	! 1-initialisation
1143	!-----------------
1144	label_max=1 ! numero de la tache, la premiere tache est notée 1
1145	label=1
1146	do i=1,n_ta_max
1147	compt_marais(i)=i
1148	enddo
1149	!$OMP PARALLEL
1150	!$OMP WORKSHARE
1151	table_out_marais(:,:) = 0
1152	marais(:) = .true.
1153	nb_pts_marais(:) = 0
1154	!$OMP END WORKSHARE
1155	!$OMP END PARALLEL
1156
1157	! 2-reperage des taches
1158	!----------------------
1159
1160	do j=2,ny-1
1161	do i=2,nx-1
1162	if ((ice(i,j).ge.1).and.flot(i,j)) then ! on est sur la glace qui flotte-------------------!
1163
1164	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
1165	! un des voisins est deja en glace
1166	mask(1) = table_out_marais(i-1,j)
1167	mask(2) = table_out_marais(i,j-1)
1168	label = label_max
1169
1170	! on determine la valeur de la tache minimun (>0) presente ds le masque
1171	do indice=1,mask_nb
1172	if (mask(indice).gt.0) label=min(label,mask(indice))
1173	enddo
1174
1175	! on fixe la valeur de la tache voisine minimun au point etudie (via label)
1176	table_out_marais(i,j)=label
1177
1178	!si un des voisins n'est pas glace alors la tache n'est pas un marais
1179	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
1180	marais(label)=.false.
1181	endif
1182
1183	! si 2 taches differentes sont dans le masque, il faut les identifier dans compt_marais
1184	! on lui affecte le numero de la tache fondamentale
1185
1186	! exemple on est sur le point X : 5 X
1187	do indice=1,mask_nb ! 20
1188	if(mask(indice).gt.label) then ! mask(2)=20 > 5
1189	compt_marais(mask(indice))=label ! compt_marais(20)=5
1190	if (.not.marais(mask(indice))) marais(label)=.false. ! si la tache n'etais pas un marais => marais =.false. marais(-(-5))=.false.
1191	where (table_out_marais(:,:).eq.mask(indice)) ! where table_out_marais(:,:)=mask(2)=20
1192	table_out_marais(:,:)=label ! table_out_marais(:,:)=label=5
1193	endwhere
1194	endif
1195	enddo
1196
1197	else !aucun des voisins est une tache
1198	table_out_marais(i,j)= label_max
1199	compt_marais(label_max)=label_max
1200
1201	! si un des voisins n'est pas glace alors la tache n'est pas un marais
1202	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
1203	marais(label_max)=.false.
1204	endif
1205	label_max = label_max+1
1206	if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches=',label_max
1207	endif
1208	else ! on est pas sur une tache--------------------------------------------
1209	table_out_marais(i,j)=0 ! Pas necessaire (reecrit 0 sur 0)
1210	endif !---------------------------------------------------------------------
1211	enddo
1212	enddo
1213
1214	marais(0)=.false.
1215
1216	!$OMP PARALLEL
1217	!$OMP DO
1218	do j=1,ny
1219	do i=1,nx
1220	if (table_out_marais(i,j).ne.0) then
1221	nb_pts_marais(compt_marais(table_out_marais(i,j)))= nb_pts_marais(compt_marais(table_out_marais(i,j)))+1
1222	endif
1223	flot_marais(i,j) = marais(table_out_marais(i,j))
1224	enddo
1225	enddo
1226	!$OMP END DO
1227	!$OMP END PARALLEL
1228
1229	debug_3D(:,:,125)=real(table_out_marais(:,:))
1230
1231	end subroutine determin_marais
1232
1233	end module flottab_mod
1234
1235

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

Download in other formats: