Context Navigation

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

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

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

Download in other formats: