Context Navigation

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

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

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

Download in other formats: