Context Navigation

flottab2-0.5-heino.f90 @ 237

Last change on this file since 237 was 237, checked in by aquiquet, 6 years ago
Sealevel is now treated as a 2D variable (sealevel_2d while sealevel remains the eustatic sea level), results should remain identical as sealevel_2d is equal to sealevel in this revision.
File size: 32.3 KB

Line
1	!> \file flottab2-0.5-heino.f90
2	!!
3	!!
4	!<
5
6	!> SUBROUTINE: massb_perturb_Tparam
7	!! Cette routine determine les endroits ou la glace
8	!! flotte , les streams et la position du front de glace
9	!! \author Vincent et Catherine
10	!! \date 10 juillet 2005
11	!! @note Il y a 4 sortes de zone
12	!! @note - Pose
13	!! @note - Grounding zone et streams gzmx et gzmy
14	!! @note - Iles ilemx, ilemy
15	!! @note - flottant flot sur le noeud majeur, flotmx sur le noeud mineur
16	!! @note Used modules:
17	!! @note - use module3D_phy
18	!! @note - use module_choix
19	!<
20
21
22
23	! =======================================================
24	!
25	SUBROUTINE FLOTTAB()
26	!
27	! Vince 5 Jan 95
28	! Modifie 20 Jan 95
29	! Modifie le 30 Novembre 98
30	! Passage f90 + determination des fronts Vincent dec 2003
31	! nettoyage et nouvelle détermination de gzmx et gzmy Cat le 10 juillet 2005
32	!
33	! -----------------
34	!
35	! Cette routine determine les endroits ou la glace
36	! flotte , les streams et la position du front de glace
37	!
38	! Il y a 4 sortes de zone
39	! Pose
40	! Grounding zone et streams gzmx et gzmy
41	! Iles ilemx, ilemy
42	! flottant flot sur le noeud majeur, flotmx sur le noeud mineur
43
44	! passage dans flottab tous les pas de temps dt :
45	!
46	! QUESTION : que se passe t'il si un noeud passe en flottant entre 2 dtt
47	!
48	! =======================================================
49
50	USE module3D_phy
51	! use printtable
52	! use sliding_heino
53	use module_choix
54
55	IMPLICIT NONE
56
57	real :: surnet !surnet hauteur de glace au dessus de la mer
58	real :: archim ! test de flottaison
59
60	logical :: gz1
61	logical :: fl1
62
63
64	real,dimension(nx,ny) :: uxs1 ! uxbar a l'entree de flottab
65	real,dimension(nx,ny) :: uys1 ! uybar a l'entree de flottab
66
67
68	integer pmx,pmy !pm=plus-moins -1 ou 1 pour x et y
69
70
71	! Variables pour la determination des differents shelfs/stream
72	! (representés comme des taches ou l'on resoud l'eq elliptique)
73	!________________________________________________________________
74	integer,parameter :: n_ta_max=2000! nombre de tache max
75	integer,dimension(0:nx,0:ny) :: table_out ! pour les numeros des taches
76	integer,dimension(n_ta_max) :: compt !contient les equivalence entre les taches
77	integer,dimension(0:n_ta_max) :: nb_pts_tache ! indique le nombre de points par tache
78	logical,dimension(0:n_ta_max) :: iceberg ! T si iceberg, F si calotte posée
79
80	! Variables pour determiner le point le plus haut (surf)
81	! d'une ile completement stream
82	!_________________________________________________________
83
84	! icetrim : T si ice stream, F si calotte posée(vertical shear)
85	logical,dimension(0:2000) :: icetrim
86
87	integer :: ii,jj
88	integer :: smax_i
89	integer :: smax_j
90	real :: smax_
91
92
93	!_________________________________________________________
94	if (itracebug.eq.1) call tracebug(' Heino: entree dans Flottab')
95
96
97
98	SHELFY = .FALSE.
99
100
101	! 1-INITIALISATION
102	! ----------------
103	! initialisation des variables pour detecter les points qui se mettent
104	! a flotter entre 2 dtt
105
106	appel_new_flot=.false.
107	do j=1,ny
108	do i=1,nx
109	new_flot_point(i,j)=.false.
110	new_flotmx(i,j)=.false.
111	new_flotmy(i,j)=.false.
112	enddo
113	enddo
114
115	! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m
116
117	ICE(:,:)=0
118	front(:,:)=0
119	frontfacex(:,:)=0
120	frontfacey(:,:)=0
121	isolx(:,:)=.FALSE.
122	isoly(:,:)=.FALSE.
123	boost=.false.
124
125	! fin de l'initialisation
126	!_____________________________________________________________________
127
128	! 2-TESTE LES NOUVEAUX POINTS FLOTTANTS
129	! -------------------------------------
130
131	do j=1,ny
132	do i=1,nx
133
134	uxs1(i,j)=uxbar(i,j)
135	uys1(i,j)=uybar(i,j)
136
137	archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j)
138
139	arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then ! le point flotte
140
141	ex_pose: if ((.not.FLOT(I,J)).and.(isynchro.eq.1)) then ! il ne flottait pas avant
142	FLOT(I,J)=.true.
143	BOOST=.false.
144
145
146	if (igrdline.eq.1) then ! en cas de grounding line prescrite
147	flot(i,j)=.false.
148	H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro ! pour avoir archim=10
149	new_flot_point(i,j)=.false.
150	endif
151
152
153
154
155	else ! isynchro=0 ou il flottait déja
156
157	if (.not.FLOT(I,J)) then ! il ne flottait pas (isynchro=0)
158	new_flot_point(i,j)=.true. ! signale un point qui ne flottait pas
159	! au pas de temps precedent
160	flot(i,j)=.true.
161
162	if (igrdline.eq.1) then ! en cas de grounding line prescrite
163	flot(i,j)=.false.
164	H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro ! pour avoir archim=10
165	new_flot_point(i,j)=.false.
166	endif
167
168	endif
169	endif ex_pose
170
171
172	else ! le point ne flotte pas
173	if(FLOT(I,J)) then ! mais il flottait avant
174	FLOT(I,J)=.false.
175	BOOST=.false.
176	endif
177	endif arch
178
179	! Si la glace flotte -> PRUDENCE !!!
180	! S et B sont alors determines avec la condition de flottabilite
181
182	if (flot(i,j)) then
183	shelfy = .true.
184
185	surnet=H(i,j)*(1.-ro/row)
186	S(i,j)=surnet+sealevel_2d(i,j)
187	B(i,j)=S(i,j)-H(i,j)
188	end if
189
190	end do
191	end do
192
193
194	IF (BOOST) GOTO 11
195	! vincent mars2004 on utilise plus boost, on redefini
196	! le front a chaque pas de temps
197	! ou alors le mettre apres la determination du front
198	! ou alors gardre la valeur des fronts : ne pas reinitialiser au
199	! debut
200
201
202
203	domain: do j=2,ny
204	do i=2,nx
205
206	! 3 A- NOUVELLE DEFINITION DE FLOTMX, LES POINTS
207	! AYANT UN DES VOISINS FLOTTANTS SONT FLOTMX ET GZMX
208	! -------------------------------------------
209
210	! fl1 est l'ancienne valeur de flotmx
211	gz1=gzmx(i,j) ! gz1 est l'ancienne valeur de gzmx
212	fl1=flotmx(i,j) ! fl1 est l'ancienne valeur de flotmx
213
214	flotmx(i,j)=flot(i,j).and.flot(i-1,j)
215
216	! test pour detecter les nouveaux flotmx entre 2 dtt :
217
218	if (flotmx(i,j).and.(new_flot_point(i,j).or. &
219	new_flot_point(i-1,j))) then
220	appel_new_flot=.true.
221	new_flotmx(i,j)=.true.
222	endif
223
224
225	! determination de gzmx
226	!__________________________________________________________________________
227
228	! gzmx si un des deux voisins est flottant et l'autre posé
229	! i-1 i
230	gzmx(i,j)=((flot(i,j).and..not.flot(i-1,j)) & ! F P
231	.or.(.not.flot(i,j).and.flot(i-1,j))) ! P F
232
233	! A condition d'etre assez proche de la flottaison
234	! sur le demi noeud condition archim > 100 m
235
236	archim=(Bsoc(i,j)+Bsoc(i-1,j))0.5-sealevel_2d(i,j)+ro/rowHmx(i,j)
237	gzmx(i,j)=gzmx(i,j).and.(archim.le.100.)
238
239	! mais gzmx peut être quand même vrai selon une conditon sur la pression effective
240	gzmx(i,j)=gzmx(i,j).or.((Neffmx(i,j).lt.neffratio*hmx(i,j)) &
241	.and..not.flotmx(i,j))
242
243
244
245	! ce cas est un ilemx traité plus loin
246	! if (i.gt.2.AND.i.lt.nx) gzmx(i,j)=gzmx(i,j).or.(.not.flot(i,j).and. &
247	! .not.flot(i-1,j).and.flot(i+1,j).and.flot(i-2,j)) ! F P P F
248
249
250	! On rajoute une condition sur la contrainte basale
251	! si tobmx(i,j) > toblim (de l'ordre de 1 bar), passage a pose
252
253	gzmx(i,j)=gzmx(i,j).and.(abs(tobmx(i,j)).lt.toblim)
254	gzmx(i,j)=gzmx(i,j).or.((hwater(i,j).gt.hwatstream).and..not.flotmx(i,j))
255
256	!-------------------------------------------------------------------------------------
257	! attention : pour expériences Heino
258	gzmx(i,j)=gzmx_heino(i,j)
259	shelfy=.true.
260	! shelfy=.false.
261	!____________________________________________________________________________________
262
263
264
265	! ATTENTION : A METTRe A LA FIN EN FAISANT PORTER SUR ilmx,flotmx gzmx
266	! si le point etait posé (non gz) et devient grounded
267	! definir la direction de la vitesse (moyenne des points)
268
269	if ((.not.gz1).and.(.not.fl1).and.gzmx(i,j).and. &
270	(i.gt.2).and.(i.lt.nx)) then
271	uxs1(i,j)=(uxbar(i+1,j)+uxbar(i-1,j))/2.
272	endif
273
274
275	! 3 B- NOUVELLE DEFINITION DE FLOTMY
276	! --------------------------------
277
278
279	gz1=gzmy(i,j) ! gz1 est l'ancienne valeur de gzmy
280	fl1=flotmy(i,j) ! fl1 est l'ancienne valeur de flotmy
281
282	flotmy(i,j)=flot(i,j).and.flot(i,j-1)
283
284	! test pour detecter les nouveaux flotmy entre 2 dtt :
285
286	if (flotmy(i,j).and.(new_flot_point(i,j).or. &
287	new_flot_point(i,j-1))) then
288	appel_new_flot=.true.
289	new_flotmy(i,j)=.true.
290	endif
291
292	! determination de gzmy
293	!__________________________________________________________________________
294
295	! gzmy si un des deux voisins est flottant et l'autre posé
296
297	gzmy(i,j)=((flot(i,j).and..not.flot(i,j-1)) &
298	.or.(.not.flot(i,j).and.flot(i,j-1)))
299
300	! A condition d'etre assez proche de la flottaison
301	! sur le demi noeud condition archim > 100 m
302
303	archim=(Bsoc(i,j)+Bsoc(i,j-1))0.5-sealevel_2d(i,j)+ro/rowHmy(i,j)
304	gzmy(i,j)=gzmy(i,j).and.(archim.le.100.)
305
306	! mais gzmx peut être quand même vrai selon une conditon sur la pression effective
307	gzmy(i,j)=gzmy(i,j).or.((Neffmy(i,j).lt.neffratio*hmy(i,j)) &
308	.and..not.flotmy(i,j))
309
310	! ce cas est un ilemy traité plus loin
311	! if (j.gt.2.AND.j.lt.ny) gzmy(i,j)=gzmy(i,j).or.(flot(i,j).and. &
312	! .not.flot(i,j-1).and.flot(i,j+1).and.flot(i,j-2))
313
314
315	! On rajoute une condition sur la contrainte basale
316	! si tobmy(i,j) > toblim (de l'ordre de 1 bar), passage a pose
317	gzmy(i,j)=gzmy(i,j).and.(abs(tobmy(i,j)).lt.toblim)
318	gzmy(i,j)=gzmy(i,j).or.((hwater(i,j).gt.hwatstream).and..not.flotmy(i,j))
319
320
321	!-------------------------------------------------------------------------------------
322	! attention : pour expériences Heino
323	gzmy(i,j)=gzmy_heino(i,j)
324	shelfy=.true.
325	! shelfy=.false.
326	!____________________________________________________________________________________
327
328
329	! ATTENTION : A METTRe A LA FIN EN FAISANT PORTER SUR ilmx,flotmx gzmx
330	! si le point etait posé (non gz) et devient grounded
331	! definir la direction de la vitesse (moyenne des points)
332
333	if (j.gt.2.AND.j.lt.ny) then
334	if ((.not.gz1).and.(.not.fl1).and.gzmy(i,j)) then
335	uys1(i,j)=(uybar(i,j+1)+uybar(i,j-1))/2.
336	endif
337
338
339	! Points du shelf proches de la grounding line (pour la fusion basale)
340	! fbm est vrai si le point est flottant mais un des voisins est pose
341	!_________________________________________________________________________
342
343	if (i.gt.2.AND.i.lt.nx) then
344	fbm(i,j)=flot(i,j).and. &
345	((.not.flot(i+1,j)).or.(.not.flot(i,j+1)) &
346	.or.(.not.flot(i-1,j)).or.(.not.flot(i,j-1)))
347	endif
348	endif
349
350	end do
351	end do domain
352
353
354	! Condition sur les bords
355
356	do i=2,nx
357	flotmx(i,1) = (flot(i,1).or.flot(i-1,1))
358	flotmy(i,1) = .false.
359	end do
360
361	do j=2,ny
362	flotmy(1,j) = (flot(1,j).or.flot(1,j-1))
363	flotmx(1,j) = .false.
364	end do
365
366	flotmx(1,1) = .false.
367	flotmy(1,1) = .false.
368
369	11 continue
370
371
372	! 4- determination des iles
373	! -------------------------
374
375	ilemx(:,:)=.false.
376	ilemy(:,:)=.false.
377
378	! selon x
379	do j=2,ny-1
380	do i=3,nx-2
381	! F G F
382	! x
383	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
384	if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. &
385	(sdx(i,j).LT.1.E-02)) then
386	ilemx(i,j)=.true.
387	ilemx(i+1,j)=.true.
388
389	! F G G F
390	! x
391	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
392	else if ((flot(i-1,j).and..not.flot(i,j) &
393	.and..not.flot(i+1,j)).and.flot(i+2,j).and. &
394	(sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then
395	ilemx(i,j)=.true.
396	ilemx(i+1,j)=.true.
397	ilemx(i+2,j)=.true.
398
399	! F G G F
400	! x
401	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
402	else if ((flot(i-2,j).and..not.flot(i-1,j) &
403	.and..not.flot(i,j)).and.flot(i+1,j).and. &
404	(sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then
405	ilemx(i-1,j)=.true.
406	ilemx(i,j)=.true.
407	ilemx(i+1,j)=.true.
408
409	! F G G G F
410	! x
411	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
412	else if ((i.lt.nx-2) &
413	.and.(flot(i-2,j).and..not.flot(i-1,j) &
414	.and..not.flot(i,j)).and..not.flot(i+1,j) &
415	.and.flot(i+2,j).and. &
416	(sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02 &
417	.and.sdx(i+1,j).LT.1.E-02)) then
418	ilemx(i-1,j)=.true.
419	ilemx(i,j)=.true.
420	ilemx(i+1,j)=.true.
421	ilemx(i+2,j)=.true.
422
423	endif
424
425	end do
426	end do
427
428	! selon y
429	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
475
476
477	! les iles n'ont pas de condition neff mais ont la condition toblim
478	do j=2,ny
479	do i=2,nx
480	ilemx(i,j)=ilemx(i,j).and.(abs(tobmx(i,j)).lt.toblim)
481	ilemy(i,j)=ilemy(i,j).and.(abs(tobmy(i,j)).lt.toblim)
482	end do
483	end do
484
485	! 6-On determine finalement la position des noeuds stream/shelf
486	! -------------------------------------------------------------
487
488	do j=2,ny-1
489	do i=2,nx-1
490
491	if (nt.ge.2) then ! pour ne pas faire ce calcul lors du premier passage
492
493	uxbar(i,j)=uxs1(i,j)
494	uybar(i,j)=uys1(i,j)
495	endif
496
497	FLGZMX(I,J)=(marine.and.(flotmx(i,j).or.gzmx(i,j).or.ilemx(i,j))) &
498	.or.(.not.marine.and.flotmx(i,j))
499	FLGZMY(I,J)=(marine.and.(flotmY(i,j).or.gzmY(i,j).or.ilemy(i,j))) &
500	.or.(.not.marine.and.flotmy(i,j))
501
502	end do
503	end do
504
505	! 5-On determine maintenant la position du front de glace
506	! -------------------------------------------------------
507	! C'est ici que l'on determine la position et le type de front
508	! print*,'on est dans flottab pour definir les fronts'
509
510
511	DO I=3,NX-2
512	DO J=3,NY-2
513	IF (H(i,j).gt.1.1) ICE(i,j)=1
514	END DO
515	END DO
516	!print*,'H ice 50,30',H(50,30),ICE(50,30)
517
518	!open(UNIT=num_file4,file='TABLE_ice-big')
519	! do i=1,nx
520	! do j=1,ny
521	! write(num_file4,*) i,j,ice(i,j)
522	! enddo
523	! enddo
524	!close (num_file4)
525	if (ISYNCHRO.eq.1) then
526	!----------------------------------------------!
527	!On determine les differents ice strean/shelf !
528	call DETERMIN_TACHE !
529	!----------------------------------------------!
530	! do i=1,nx
531	! do j=1,ny
532	! if (ice(i,j)==1.and.table_out(i,j)==0) then
533	! print*,'i,j AAAAAAAAAAAAAAA',i,j
534	! endif
535	! enddo
536	! enddo
537	ICE=0
538	!On compte comme englacé uniquement les calottes dont une partie et posé
539	! pause
540	! DO I=2,NX-1
541	! DO J=2,NY-1
542	DO I=3,NX-2
543	DO J=3,NY-2
544	IF (.not.iceberg(table_out(i,j))) THEN !!! on est pas sur un iceberg TEST1
545	if (nb_pts_tache(table_out(i,j)).ge.1) then
546	ICE(i,j)=1
547	if (nb_pts_tache(table_out(i,j)).le.5) then
548	FLGZMX(I,J)=.false.
549	FLGZMY(I,J)=.false.
550	! print*,'i,jcoupés',FLGZMX(I,J),FLGZMY(I,J),i,j
551	endif
552
553	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ici on est sur une tache non iceberg >= 5 points
554	!!!!!!!!! on teste si la tache n'est pas completement ice stream
555	if (icetrim(table_out(i,j))) THEN ! on a une ile d'ice stream !!!! TEST2
556	smax_i=0 ; smax_j=0 ; smax_=sealevel !afq -- warning: what about local sealevel?
557	DO II=3,NX-2
558	DO JJ=3,NY-2
559	if (table_out(ii,jj).eq.table_out(i,j)) then
560	if (S(ii,jj).gt.smax_) then
561	smax_ =S(ii,jj)
562	smax_i=ii
563	smax_j=jj
564	endif
565	endif
566	END DO
567	END DO
568	gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false.
569	gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false.
570	flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false.
571	flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false.
572	endif !!! !!!!!!! TEST2
573	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
574	ELSE !!! on est sur un iceberg !!!!!!!! TEST1
575	ICE(i,j)=0
576	! print*,'iceberg',i,j,h(i,j)
577	H(i,j)=1.
578	ENDIF !!!!!!!!!!!! !!!!!!!! TEST1
579	endif
580	! IF (.not.iceberg(table_out(i,j)).and.nb_pts_tache(table_out(i,j)).ge.1) THEN
581	! ICE(i,j)=1
582	! if (nb_pts_tache(table_out(i,j)).le.5) then
583	! FLGZMX(I,J)=.false.
584	! FLGZMY(I,J)=.false.
585	!print*,'i,j coupés',FLGZMX(I,J),FLGZMY(I,J),i,j
586	! endif
587	! ELSE
588	!IF (table_out(i,j).gt.0) print*,'i,j coupés',i,j,table_out(i,j),nb_pts_tache(table_out(i,j)),flgzmx(i,j),flotmx(i,j)
589	! ENDIF
590	END DO
591	END DO
592
593	!call printtable_i(ice,nom_table)
594	!ause
595	!----------------------------------------------!
596	!On caracterise le front des ice shelfs/streams!
597	call DETERMIN_FRONT !
598	!----------------------------------------------!
599
600
601	!nom_table='ice_af'
602	!call printtable_i(ice,nom_table)
603	!nom_table='GZMX2'
604	!call printtable_l(FLGZMX,nom_table)
605	!nom_table='flot'
606	!call printtable_l(FLOT,nom_table)
607	!nom_table='mk'
608	!call printtable_i(mk,nom_table)
609	!nom_table='mk0'
610	!call printtable_i(mk0,nom_table)
611	endif
612
613
614	!fin de routine flottab2
615	!print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j)
616	contains
617	!====================================================================
618	!====================================================================
619	SUBROUTINE DETERMIN_TACHE
620	implicit none
621	integer :: indice
622	integer :: label ! no des taches rencontrées dans le mask
623	integer :: label_max ! no temporaire maxi de tache rencontrées
624	! integer :: mask_nb = 4
625	integer :: mask_nb = 2 ! version ou on ne compte pas les diagonales
626	integer :: vartemp ! variable temporaire pour reorganiser compt
627	! integer,dimension(mask_nb) :: mask
628	integer,dimension(2) :: mask
629
630	! 1-initialisation
631	!-----------------
632	label_max=1 ! numero de la tache, la premiere tache est notée 1
633	label=1
634	do i=1,n_ta_max
635	compt(i)=i
636	enddo
637	! table_in = .false.
638	table_out = 0
639	iceberg = .true.
640	icetrim = .true.
641	nb_pts_tache = 0
642	! open(unit=100,file="tache.data",status='replace')
643
644	! 2-reperage des taches
645	!----------------------
646	do i=1,nx
647	do j=1,ny
648
649	if (label_max==0) print*,'labmax=0',i,j
650	IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------!
651	! if ((ice(i-1,j-1).ge.1).or.(ice(i-1,j).ge.1).or.& !masque de 4 cases adjacentes
652	! (ice(i-1,j+1).ge.1).or.(ice(i,j-1).ge.1)) then
653	! mask(1) = table_out(i-1,j-1)
654	! mask(2) = table_out(i-1,j)
655	! mask(3) = table_out(i-1,j+1)
656	! mask(4) = table_out(i,j-1)
657	if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then !masque de 2 cases adjacentes
658	! un des voisins est deja en glace
659	mask(1) = table_out(i-1,j)
660	mask(2) = table_out(i,j-1)
661	label = label_max
662
663	!on determine la valeur de la tache minimun (>0) presente ds le masque
664	do indice=1,mask_nb
665	if (mask(indice).gt.0) label=min(label,mask(indice))
666	enddo
667
668	!on fixe la valeur de la tache voisine minimun au point etudié(via label)
669	table_out(i,j)=label
670	!si ce noeud est posé, alors le tache n'est pas un iceberg et iceberg=.F.
671	if (.not.FLOT(I,J)) then
672	iceberg(label)=.false.
673	endif
674	!si ce noeud est posé, alors la tache n'est pas un ice stream et icestrim=.F.
675	if (.not.gzmx(I,J).AND..not.gzmy(I,J)) then
676	icetrim(label)=.false.
677	endif
678
679	!si 2 taches differentes sont dans le masque, il faut les identifier dans compt
680	!i.e. les plus grands numeros correspondent au plus petit
681	!on lui affecte le numero de la tache fondamentale avec un signe - pour indiquer le changement
682	do indice=1,mask_nb
683	if(mask(indice).gt.label) then
684	compt(mask(indice))=-label
685	endif
686	enddo
687
688	else !aucun des voisins est une tache
689	table_out(i,j)= label_max
690	compt(label_max)=label_max
691	!si ce noeud est posé, alors la ache n'est pas un iceberg et iceberg=.F.
692	if (.not.FLOT(I,J)) then
693	iceberg(label_max)=.false.
694	endif
695	!si ce noeud est posé, alors le tache n'est pas un ice stream et icestrim=.F.
696	if (.not.gzmx(I,J).AND..not.gzmy(I,J)) then
697	icetrim(label)=.false.
698	endif
699
700	label_max = label_max+1
701	if (label_max.gt.n_ta_max) print*,'trop de taches=',label_max
702	endif
703
704	ELSE !on est pas sur une tache----------------------------------------------!
705	table_out(i,j)=0 ! Pas necessaire (reecrit 0 sur 0) !
706	ENDIF!----------------------------------------------------------------------!
707
708
709	enddo
710	enddo
711
712	! On reorganise compt en ecrivant le numero de la tache fondamentale
713	! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité)
714	! On indique aussi le nb de point que contient chaque taches (nb_pts_tache)
715	do indice=1,label_max
716	vartemp = compt(indice)
717	if (compt(indice).lt.0) then
718	compt(indice)= compt(-vartemp)
719	if (.not.iceberg(indice)) iceberg(-vartemp)=.false.
720	if (.not.icetrim(indice)) icetrim(-vartemp)=.false.
721	endif
722
723	enddo
724	do i=1,nx
725	do j=1,ny
726	if (table_out(i,j).ne.0) then
727	table_out(i,j)=compt(table_out(i,j))
728	nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1
729	endif
730	enddo
731	enddo
732	! 3-visualisation
733	!----------------
734	!print*,'_pts_tache0:',nb_pts_tache(0:label_max+15)
735	!print*,'lab tache', label_max
736	!print*,'compt',compt
737	!om_table='tache'
738	!call printtable_i(table_out,nom_table)
739	!nom_table='bergs'
740	!call printtable_l(iceberg(table_out),nom_table)
741	!nom_table='pt_ta'
742	!call printtable_i(nb_pts_tache(table_out),nom_table)
743
744
745	END SUBROUTINE DETERMIN_TACHE
746	!====================================================================
747	!====================================================================
748	SUBROUTINE DETERMIN_FRONT
749
750	DO I=3,NX-2
751	DO J=3,NY-2
752
753	surice:IF (ICE(i,j)==0) THEN
754	do pmx=-1,1,2
755	do pmy=-1,1,2
756	diagice : if (ICE(i+pmx,j+pmy)==1) then
757	!!!!! print*,'========',i,j,ICE(i+pmx,j)+ICE(i,j+pmy)
758	if ((ICE(i+pmx,j)+ICE(i,j+pmy)==2)) then
759	! test (i) pour eviter les langues de glaces diagonales en coin(26dec04)
760	if ((ICE(i+2pmx,j)==1.and.ICE(i+2pmx,j+pmy)==0).or.&
761	(ICE(i,j+2pmy)==1.and.ICE(i+pmx,j+2pmy)==0)) then
762	! print*,i,j,'langues en coin'
763	ICE(i,j)=1
764	H(i,j)=max(1.,H(i,j))
765	endif
766	! test (i) pour eviter les langues de glaces diagonales
767	! mouvement du cheval aux echecs
768	if ((ICE(i+2pmx,j+pmy)+ICE(i+pmx,j+2pmy)==1)) then
769	if (ICE(i+2*pmx,j+pmy)==1.and. &
770	(ICE(i+2pmx,j+2pmy)+ICE(i,j+2*pmy)).ge.1) then
771	ICE(i,j)=1
772	H(i,j)=max(1.,H(i,j))
773	endif
774	if (ICE(i+pmx,j+2*pmy)==1.and. &
775	(ICE(i+2pmx,j+2pmy)+ICE(i+2*pmx,j)).ge.1) then
776	ICE(i,j)=1
777	H(i,j)=max(1.,H(i,j))
778	endif
779	! test (ii) pour eviter les langues de glaces diagonales
780	! le point glace ICE(i+pmx,j+pmy) a ses 4 voisins frontaux en glace
781	! mais 2 voisins vides diagonalement opposes
782	elseif ((ICE(i+2pmx,j+pmy)+ICE(i+pmx,j+2pmy)==2) &
783	.and.ICE(i+2pmx,j+2pmy)==0) then
784	! test (iii) pour faire les tests (i) et (ii)
785	! if ((ICE(i+2pmx,j+pmy)+ICE(i+pmx,j+2pmy).ge.1) &
786	! .and.ICE(i+2pmx,j+2pmy)==0) then
787	ICE(i,j)=1
788	H(i,j)=max(1.,H(i,j))
789	ICE(i+2pmx,j+2pmy)=1
790	H(i+2pmx,j+2pmy)=max(1.,H(i+2pmx,j+2pmy))
791	endif
792	endif
793	endif diagice
794	enddo
795	enddo
796	ENDIF surice
797	END DO
798	END DO
799	! print*,'dans remplissage baies',time
800	DO I=3,NX-2
801	DO J=3,NY-2
802	surice2:IF (ICE(i,j)==0) THEN
803	! test (iii) pour trouver les baies de largeur 1 ou 2 cases
804	! et combler les trous si ce sont des baies
805	! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles
806	! baies horizontales
807	if (ICE(i-1,j)==1.and.(ICE(i+1,j)==1.or.ICE(i+2,j)==1)) then
808	if (ICE(i,j-1)==1.or.ICE(i,j+1)==1) then! ICE(i,j)=1
809	ICE(i,j)=1
810	H(i,j)=max(1.,H(i,j))
811	endif
812	endif
813	! baies verticales
814	if (ICE(i,j-1)==1.and.(ICE(i,j+1)==1.or.ICE(i+2,j)==1)) then
815	if (ICE(i-1,j)==1.or.ICE(i+1,j)==1) then !ICE(i,j)=1
816	ICE(i,j)=1
817	H(i,j)=max(1.,H(i,j))
818	endif
819	! ICE(i,j)=1
820	endif
821	! if (ICE(i-1,j)==1.and.(ICE(i+1,j)==1.or.ICE(i+2,j)==1)) ICE(i,j)=1
822	! if (ICE(i,j-1)==1.and.(ICE(i,j+1)==1.or.ICE(i+2,j)==1)) ICE(i,j)=1
823	ENDIF surice2
824	END DO
825	END DO
826	DO I=3,NX-2
827	DO J=3,NY-2
828	surice3:IF (ICE(i,j)==0) THEN
829	!! test (iii) pour trouver les baies de largeur 1 ou 2 cases (oublie au test d'avant)
830	!! et combler les trous si ce sont des baies
831	!! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles
832	!! baies horizontales
833	if (ICE(i-1,j)==1.and.(ICE(i+1,j)==1.or.ICE(i+2,j)==1)) then
834	if (ICE(i,j-1)==1.or.ICE(i,j+1)==1) then!ICE(i,j)=1
835	ICE(i,j)=1
836	H(i,j)=max(1.,H(i,j))
837	! if (ISYNCHRO.eq.1) print*,'surice3 hor',i,j
838	endif
839	! ICE(i,j)=1
840	endif
841	!! baies verticales
842	if (ICE(i,j-1)==1.and.(ICE(i,j+1)==1.or.ICE(i+2,j)==1)) then
843	if (ICE(i-1,j)==1.or.ICE(i+1,j)==1) then!ICE(i,j)=1
844	ICE(i,j)=1
845	H(i,j)=max(1.,H(i,j))
846	! if (ISYNCHRO.eq.1) print*,'surice3 ver',i,j
847	endif
848
849	! ICE(i,j)=1
850	endif
851	!! if (ICE(i-1,j)==1.and.(ICE(i+1,j)==1.or.ICE(i+2,j)==1)) ICE(i,j)=1
852	!! if (ICE(i,j-1)==1.and.(ICE(i,j+1)==1.or.ICE(i+2,j)==1)) ICE(i,j)=1
853	ENDIF surice3
854	END DO
855	END DO
856
857	DO I=1,NX
858	DO J=1,NY
859
860	IF (ICE(i,j)==1) THEN !!!!!!!!!!!!!!!!!!test SI ICE=1
861	! if ICE, on determine FRONT...
862	! ainsi, FRONT=0 sur les zones = 0
863	FRONT(i,j)=(ICE(i-1,j)+ICE(i+1,j)+ICE(i,j+1)+ICE(i,j-1))
864	!front= le nb de faces en contact avec un voisin englacé
865	ENDIF
866	END DO
867	END DO
868
869	DO I=1,NX
870	DO J=1,NY
871	!On ne compte pas les icebergs de 2 cases (horizontales ou verticales)
872	IF (FRONT(i,j).eq.1) then
873	IF (FRONT(i+1,j).eq.1) then
874	ICE(i,j)=0
875	ICE(i+1,j)=0
876	FRONT(i,j)=0
877	FRONT(i+1,j)=0
878	ENDIF
879	IF (FRONT(i,j+1).eq.1) then
880	ICE(i,j)=0
881	ICE(i,j+1)=0
882	FRONT(i,j)=0
883	FRONT(i,j+1)=0
884	ENDIF
885	ENDIF
886
887	IF (FRONT(i,j).ge.1.and.FRONT(i,j).le.3) THEN !front(entre 1 et 3)
888
889	if ((ICE(i-1,j)+ICE(i+1,j)).lt.2) then
890	! il y a un front // a x
891	if ((ICE(i-1,j)+ICE(i+1,j)).eq.0) then
892	isolx(i,j)=.TRUE.
893	elseif (ICE(i-1,j).eq.0) then
894	frontfacex(i,j)=-1 ! front i-1 \|i i+1
895	else
896	frontfacex(i,j)=+1 ! front i-1 i\| i+1
897	endif
898	endif
899
900
901	if ((ICE(i,j-1)+ICE(i,j+1)).lt.2) then
902	! il y a un front // a y
903	if ((ICE(i,j-1)+ICE(i,j+1)).eq.0) then
904	isoly(i,j)=.TRUE. !front j-1 \|j\| j+1
905	elseif (ICE(i,j-1).eq.0) then
906	frontfacey(i,j)=-1 !front j-1 \|j j+1
907	else
908	frontfacey(i,j)=+1 !front j-1 j\| j+1
909	endif
910	endif
911	!isolx signifie pas de voisins en x
912	!isoly signifie pas de voisins en y
913	!Remarque :
914	!si isolx/y=.true. alors frontfacex/y=0 (a la fois +1 & -1 or +1-1=0)
915
916	END IF !fin du test il y a un front
917	!!!!!!!!!!!!!!!! END IF !!!!!!!!!!!!!!!!!!fin du test SI ICE=1
918	END DO
919	END DO
920	END SUBROUTINE DETERMIN_FRONT
921	!====================================================================
922	!====================================================================
923
924	END SUBROUTINE FLOTTAB
925

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

Context Navigation

source: trunk/SOURCES/Heino_files/flottab2-0.5-heino.f90 @ 237

Download in other formats: