Context Navigation

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

Last change on this file since 65 was 58, checked in by aquiquet, 8 years ago
Bug fix in flottab: surface below sea level after isostasy
File size: 34.6 KB

Line
1	!> \file flottab2-0.7.f90
2	!! Module pour determiner les endroits ou la glace
3	!! flotte , les iles, et la position du front de glace
4	!<
5
6	!> \namespace flottab_mod
7	!! Determine les endroits ou la glace
8	!! flotte , les iles, et la position du front de glace
9	!! \author Vincent & Cat
10	!! \date 10 juillet 2005
11	!! @note Used module
12	!! @note - use module3D_phy
13	!! @note - use module_choix
14	!! @todo "nasty Island". If the bedrock is above sealevel force grounded (mk_init)
15	!<
16	module flottab_mod
17
18	USE module3D_phy
19	use module_choix
20
21
22
23	IMPLICIT NONE
24
25	real :: surnet !< surnet hauteur de glace au dessus de la mer
26	real :: archim !< test de flottaison
27
28	integer:: itestf
29
30	logical,dimension(nx,ny) :: gz1mx,gz1my
31	logical,dimension(nx,ny) :: fl1mx,fl1my
32
33
34	real,dimension(nx,ny) :: uxs1 !< uxbar a l'entree de flottab
35	real,dimension(nx,ny) :: uys1 !< uybar a l'entree de flottab
36
37
38	integer pmx,pmy !pm=plus-moins -1 ou 1 pour x et y
39
40
41	! Variables pour la determination des differents shelfs/stream
42	! (representés comme des taches ou l'on resoud l'eq elliptique)
43	!________________________________________________________________
44	integer,parameter :: n_ta_max=2000!< nombre de tache max
45	integer,dimension(nx,ny) :: table_out !< pour les numeros des taches
46	integer,dimension(nx,ny) :: tablebis !< pour les numeros des taches
47	integer,dimension(0:n_ta_max) :: compt !< contient les equivalence entre les taches
48	integer,dimension(0:n_ta_max) :: nb_pts_tache !< indique le nombre de points par tache
49	logical,dimension(0:n_ta_max) :: iceberg !< T si iceberg, F si calotte posee
50
51	logical,dimension(nx,ny) :: mask_tache_ij !< masque de travail
52	!< vrai pour toute la tache de i,j
53	integer,dimension(2) :: smax_coord !< pour le maxloc des iles
54
55	! Variables pour determiner le point le plus haut (surf)
56	! d'une ile completement stream
57	!_________________________________________________________
58
59	! icetrim : T si ice stream, F si calotte posee(vertical shear)
60	logical,dimension(n_ta_max) :: icetrim
61
62	integer :: ii,jj
63	integer :: smax_i
64	integer :: smax_j
65	real :: smax_
66	integer :: numtache
67	integer :: nb_pt
68	contains
69	! -----------------------------------------------------------------------------------
70	!> SUBROUTINE: flottab()
71	!! Cette routine determine les endroits ou la glace
72	!! flotte , les iles, et la position du front de glace
73	!! @note Il y a 4 sortes de zone
74	!! @note - Pose
75	! @note - Grounding zone et streams gzmx et gzmy
76	! @note - Iles ilemx, ilemy
77	! @note - flottant flot sur le noeud majeur, flotmx sur le noeud mineur
78	!>
79	subroutine flottab()
80	!
81	! Vince 5 Jan 95
82	! Modifie 20 Jan 95
83	! Modifie le 30 Novembre 98
84	! Passage f90 + determination des fronts Vincent dec 2003
85	! nettoyage et nouvelle détermination de gzmx et gzmy Cat le 10 juillet 2005
86	! Re-nettoyage par Cat en aout 2006.
87	! Le calcul de gzmx et gzmy pour les points intérieurs passe dans la subroutine dragging
88	! pour les points cotiers, toujours fait dans flottab
89	!
90	! -----------------
91	!
92	! Cette routine determine les endroits ou la glace
93	! flotte , les iles, et la position du front de glace
94	!
95	! Il y a 4 sortes de zone
96	! Pose
97	! Grounding zone et streams gzmx et gzmy
98	! Iles ilemx, ilemy
99	! flottant flot sur le noeud majeur, flotmx sur le noeud mineur
100
101	! passage dans flottab tous les pas de temps dt !
102	!
103	! _________________________________________________________
104
105
106
107	if (itracebug.eq.1) call tracebug(' Entree dans routine flottab')
108
109	SHELFY = .FALSE.
110
111
112	! cas particulier des runs paleo ou on impose un masque grounded
113
114	if (igrdline.eq.2) then
115	where ( mk_init(:,:).eq.1) ! pose
116	flot(:,:) = .False.
117	H(:,:)=max(H(:,:),(10.+sealevel-Bsoc(:,:))*row/ro) ! pour avoir archim=10
118	S(:,:) = Bsoc(:,:) + H(:,:)
119	elsewhere
120	flot(:,:) = .True.
121	end where
122	end if
123
124
125	! 1-INITIALISATION
126	! ----------------
127	! initialisation des variables pour detecter les points qui se mettent
128	! a flotter entre 2 dtt
129
130	appel_new_flot=.false.
131	do j=1,ny
132	do i=1,nx
133	new_flot_point(i,j)=.false.
134	new_flotmx(i,j)=.false.
135	new_flotmy(i,j)=.false.
136	enddo
137	enddo
138
139	! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m
140
141	ICE(:,:)=0
142	front(:,:)=0
143	frontfacex(:,:)=0
144	frontfacey(:,:)=0
145	isolx(:,:)=.FALSE.
146	isoly(:,:)=.FALSE.
147	cotemx(:,:)=.false.
148	cotemy(:,:)=.false.
149	boost=.false.
150
151	! fin de l'initialisation
152	!_____________________________________________________________________
153
154	! 2-TESTE LES NOUVEAUX POINTS FLOTTANTS
155	! -------------------------------------
156
157
158	do j=1,ny
159	do i=1,nx
160
161	uxs1(i,j)=uxbar(i,j)
162	uys1(i,j)=uybar(i,j)
163
164	archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
165
166
167	arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then ! le point flotte
168	mk(i,j)=1
169
170
171	ex_pose: if ((.not.FLOT(I,J)).and.(isynchro.eq.1)) then ! il ne flottait pas avant
172	FLOT(I,J)=.true.
173	BOOST=.false.
174
175	! Attention : avec le bloc dessous il faut faire un calcul de flottaison a la lecture
176
177	if (igrdline.eq.1) then ! en cas de grounding line prescrite
178	flot(i,j)=.false.
179	H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro ! pour avoir archim=10
180	new_flot_point(i,j)=.false.
181	endif
182
183
184
185
186	else ! isynchro=0 ou il flottait déja
187
188	if (.not.FLOT(I,J)) then ! il ne flottait pas (isynchro=0)
189	new_flot_point(i,j)=.true. ! signale un point qui ne flottait pas
190	! au pas de temps precedent
191	flot(i,j)=.true.
192
193	if (igrdline.eq.1) then ! en cas de grounding line prescrite
194	flot(i,j)=.false.
195	H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro ! pour avoir archim=10
196	new_flot_point(i,j)=.false.
197	endif
198
199	endif
200	endif ex_pose
201
202
203	else ! le point ne flotte pas
204	mk(i,j)=0
205
206	if(FLOT(I,J)) then ! mais il flottait avant
207	FLOT(I,J)=.false.
208	BOOST=.false.
209	endif
210	endif arch
211
212	! Si la glace flotte -> PRUDENCE !!!
213	! S et B sont alors determines avec la condition de flottabilite
214
215	if (flot(i,j)) then
216	shelfy = .true.
217
218	surnet=H(i,j)*(1.-ro/row)
219	S(i,j)=surnet+sealevel
220	B(i,j)=S(i,j)-H(i,j)
221	end if
222
223	end do
224	end do
225
226
227	!!$ do i=1,nx
228	!!$ do j=1,ny
229	!!$ if (flot(i,j)) then
230	!!$ mk(i,j)=1
231	!!$ else
232	!!$ mk(i,j)=0
233	!!$ endif
234	!!$ end do
235	!!$ end do
236
237
238
239
240
241	!-----------------------------------------------------------------------
242
243	domain_x: do j=1,ny
244	do i=2,nx
245
246	! 3_x A- NOUVELLE DEFINITION DE FLOTMX, LES POINTS
247	! AYANT UN DES VOISINS FLOTTANTS SONT FLOTMX ET GZMX
248	! -------------------------------------------
249
250	! fl1 est l'ancienne valeur de flotmx
251	gz1mx(i,j)=gzmx(i,j) ! gz1 est l'ancienne valeur de gzmx
252	fl1mx(i,j)=flotmx(i,j) ! fl1 est l'ancienne valeur de flotmx
253
254	flotmx(i,j)=flot(i,j).and.flot(i-1,j)
255
256	! test pour detecter les nouveaux flotmx entre 2 dtt :
257
258	if (flotmx(i,j).and.(new_flot_point(i,j).or. &
259	new_flot_point(i-1,j))) then
260	appel_new_flot=.true.
261	new_flotmx(i,j)=.true.
262	endif
263
264
265	! premiere determination de gzmx
266	!__________________________________________________________________________
267
268	! gzmx si un des deux voisins est flottant et l'autre posé
269	! i-1 i
270	gzmx(i,j)=((flot(i,j).and..not.flot(i-1,j)) & ! F P
271	.or.(.not.flot(i,j).and.flot(i-1,j))) ! P F
272
273	! A condition d'etre assez proche de la flottaison
274	! sur le demi noeud condition archim < 100 m
275
276	archim=(Bsoc(i,j)+Bsoc(i-1,j))0.5-sealevel+ro/rowHmx(i,j)
277	gzmx(i,j)=gzmx(i,j).and.(archim.le.100.)
278	cotemx(i,j)=gzmx(i,j)
279
280	end do
281	end do domain_x
282
283	if (itracebug.eq.1) call tracebug(' routine flottab apres domain_x')
284
285	! 3_y B- NOUVELLE DEFINITION DE FLOTMY
286	! --------------------------------
287	domain_y: do j=2,ny
288	do i=1,nx
289
290	gz1my(i,j)=gzmy(i,j) ! gz1 est l'ancienne valeur de gzmy
291	fl1my(i,j)=flotmy(i,j) ! fl1 est l'ancienne valeur de flotmy
292
293	flotmy(i,j)=flot(i,j).and.flot(i,j-1)
294
295	! test pour detecter les nouveaux flotmy entre 2 dtt :
296
297	if (flotmy(i,j).and.(new_flot_point(i,j).or. &
298	new_flot_point(i,j-1))) then
299	appel_new_flot=.true.
300	new_flotmy(i,j)=.true.
301	endif
302
303	! premiere determination de gzmy
304	!__________________________________________________________________________
305
306	! gzmy si un des deux voisins est flottant et l'autre posé
307
308	gzmy(i,j)=((flot(i,j).and..not.flot(i,j-1)) &
309	.or.(.not.flot(i,j).and.flot(i,j-1)))
310
311	! A condition d'etre assez proche de la flottaison
312	! sur le demi noeud condition archim > 100 m
313
314	archim=(Bsoc(i,j)+Bsoc(i,j-1))0.5-sealevel+ro/rowHmy(i,j)
315	gzmy(i,j)=gzmy(i,j).and.(archim.le.100.)
316	cotemy(i,j)=gzmy(i,j)
317
318	end do
319	end do domain_y
320
321
322
323	!-------------------------------------------------------------------------------------
324	! attention : pour expériences Heino
325	! gzmy(i,j)=gzmy_heino(i,j)
326	! shelfy=.true.
327	! shelfy=.false.
328	!____________________________________________________________________________________
329
330
331	!!$
332	!!$
333	!!$! 4- Condition sur les bords
334	!!$
335	!!$
336	!!$ do i=2,nx
337	!!$ flotmx(i,1) = (flot(i,1).or.flot(i-1,1))
338	!!$ flotmy(i,1) = .false.
339	!!$ end do
340	!!$
341	!!$ do j=2,ny
342	!!$ flotmy(1,j) = (flot(1,j).or.flot(1,j-1))
343	!!$ flotmx(1,j) = .false.
344	!!$ end do
345	!!$
346	!!$ flotmx(1,1) = .false.
347	!!$ flotmy(1,1) = .false.
348	!!$
349
350
351	! 4- determination des iles
352	! -------------------------
353
354	ilemx(:,:)=.false.
355	ilemy(:,:)=.false.
356
357	! selon x
358	ilesx: do j=2,ny-1
359	do i=3,nx-2
360	! F G F
361	! x
362	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
363	if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. &
364	(sdx(i,j).LT.1.E-02)) then
365	ilemx(i,j)=.true.
366	ilemx(i+1,j)=.true.
367
368	! F G G F
369	! x
370	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
371	else if ((flot(i-1,j).and..not.flot(i,j) &
372	.and..not.flot(i+1,j)).and.flot(i+2,j).and. &
373	(sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then
374	ilemx(i,j)=.true.
375	ilemx(i+1,j)=.true.
376	ilemx(i+2,j)=.true.
377
378	! F G G F
379	! x
380	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
381	else if ((flot(i-2,j).and..not.flot(i-1,j) &
382	.and..not.flot(i,j)).and.flot(i+1,j).and. &
383	(sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then
384	ilemx(i-1,j)=.true.
385	ilemx(i,j)=.true.
386	ilemx(i+1,j)=.true.
387
388	! F G G G F
389	! x
390	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
391	else if ((i.lt.nx-2) &
392	.and.(flot(i-2,j).and..not.flot(i-1,j) &
393	.and..not.flot(i,j)).and..not.flot(i+1,j) &
394	.and.flot(i+2,j).and. &
395	(sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02 &
396	.and.sdx(i+1,j).LT.1.E-02)) then
397	ilemx(i-1,j)=.true.
398	ilemx(i,j)=.true.
399	ilemx(i+1,j)=.true.
400	ilemx(i+2,j)=.true.
401
402	endif
403
404	end do
405	end do ilesx
406
407	! selon y
408	ilesy: do j=3,ny-2
409	do i=2,nx-1
410	! F G F
411	! x
412	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
413	if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. &
414	(sdy(i,j).LT.1.E-02)) then
415	ilemy(i,j)=.true.
416	ilemy(i,j+1)=.true.
417
418	! F G G F
419	! x
420	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
421	else if ((flot(i,j-1).and..not.flot(i,j) &
422	.and..not.flot(i,j+1)).and.flot(i,j+2).and. &
423	(sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then
424	ilemy(i,j)=.true.
425	ilemy(i,j+1)=.true.
426	ilemy(i,j+2)=.true.
427
428	! F G G F
429	! x
430	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
431	else if ((flot(i,j-2).and..not.flot(i,j-1) &
432	.and..not.flot(i,j)).and.flot(i,j+1).and. &
433	(sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then
434	ilemy(i,j-1)=.true.
435	ilemy(i,j)=.true.
436	ilemy(i,j+1)=.true.
437
438	! F G G G F
439	! x
440	! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
441	else if ((j.lt.ny-2) &
442	.and.(flot(i,j-2).and..not.flot(i,j-1) &
443	.and..not.flot(i,j)).and..not.flot(i,j+1) &
444	.and.flot(i,j+2).and. &
445	(sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02 &
446	.and.sdy(i,j+1).LT.1.E-02)) then
447	ilemy(i,j-1)=.true.
448	ilemy(i,j)=.true.
449	ilemy(i,j+1)=.true.
450	ilemy(i,j+2)=.true.
451	endif
452	end do
453	end do ilesy
454	! fin des iles
455
456	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
457	!!$if (itestf.gt.0) then
458	!!$ write(6,*) 'dans flottab avant dragging asymetrie sur H pour time=',time
459	!!$ stop
460	!!$else
461	!!$ write(6,*) 'dans flottab avant dragging pas d asymetrie sur H pour time=',time
462	!!$
463	!!$end if
464
465	! 5- calcule les noeuds qui sont streams a l'interieur et donne le betamx et betamy
466	!----------------------------------------------------------------------------------
467	if (iter_beta.eq.0) then
468	Call dragging
469	endif
470	if (itracebug.eq.1) call tracebug(' routine flottab apres call dragging')
471	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
472	!!$if (itestf.gt.0) then
473	!!$ write(6,*) 'dans flottab apres dragging asymetrie sur H pour time=',time
474	!!$ stop
475	!!$else
476	!!$ write(6,*) 'dans flottab aapres dragging pas d asymetrie sur H pour time=',time
477	!!$
478	!!$end if
479
480	! 6- calcule les vitesses des points qui sont devenus gzm
481
482	do j=1,ny
483	do i=2,nx-1
484	! si le point etait posé (non gz) et devient gzmx
485	! definir la direction de la vitesse (moyenne des points)
486
487	if ((.not.gz1mx(i,j)).and.(.not.fl1mx(i,j)).and.gzmx(i,j).and. &
488	(i.gt.2).and.(i.lt.nx)) then
489	uxs1(i,j)=(uxbar(i+1,j)+uxbar(i-1,j))/2.
490	endif
491
492	end do
493	end do
494
495	do j=2,ny-1
496	do i=1,nx
497	! si le point etait posé (non gz) et devient gzmy
498	! definir la direction de la vitesse (moyenne des points)
499	if ((.not.gz1my(i,j)).and.(.not.fl1my(i,j)).and.gzmy(i,j).and. &
500	(j.gt.2).and.(j.lt.ny)) then
501	uys1(i,j)=(uybar(i,j+1)+uybar(i,j-1))/2.
502	endif
503
504	end do
505	end do
506
507
508
509	! 7-On determine finalement la position des noeuds stream ou shelf
510	! -------------------------------------------------------------
511
512	if (nt.ge.2) then ! pour ne pas faire ce calcul lors du premier passage
513	uxbar(:,:)=uxs1(:,:)
514	uybar(:,:)=uys1(:,:)
515	endif
516
517	flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:))) &
518	.or.(.not.marine.and.flotmx(:,:))
519	flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:))) &
520	.or.(.not.marine.and.flotmy(:,:))
521
522
523
524	! 8- Pour la fusion basale sous les ice shelves- region proche de la grounding line
525	!---------------------------------------------------------------------------------
526	! fbm est vrai si le point est flottant mais un des voisins est pose
527	!_________________________________________________________________________
528	do j=2,ny-1
529	do i=2,nx-1
530
531	! if (i.gt.2.AND.i.lt.nx) then
532	fbm(i,j)=flot(i,j).and. &
533	((.not.flot(i+1,j)).or.(.not.flot(i,j+1)) &
534	.or.(.not.flot(i-1,j)).or.(.not.flot(i,j-1)))
535	! endif
536	end do
537	end do
538
539
540
541	! 9-On determine maintenant la position du front de glace
542	! -------------------------------------------------------
543	! C'est ici que l'on determine la position et le type de front
544	! print*,'on est dans flottab pour definir les fronts'
545
546
547
548	!!$do i=3,nx-2
549	!!$ do j=3,ny-2
550	!!$ if (h(i,j).gt.1.1) ice(i,j)=1
551	!!$ end do
552	!!$end do
553
554	where (flot(:,:))
555	where (H(:,:).gt.(1.1))
556	ice(:,:)=1
557	elsewhere
558	ice(:,:)=0
559	end where
560	elsewhere
561	where (H(:,:).gt.(.1))
562	ice(:,:)=1
563	elsewhere
564	ice(:,:)=0
565	end where
566	end where
567
568	call DETERMIN_TACHE
569
570	synchro : if (isynchro.eq.1) then
571	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
572	!!$if (itestf.gt.0) then
573	!!$ write(6,*) 'dans flottab avant DETERMIN_TACHE asymetrie sur H pour time=',time
574	!!$ stop
575	!!$else
576	!!$ write(6,*) 'dans flottab apres DETERMIN_TACHE pas d asymetrie sur H pour time=',time
577	!!$
578	!!$end if
579
580
581	!----------------------------------------------!
582	!On determine les differents ice strean/shelf !
583	call DETERMIN_TACHE !
584	!----------------------------------------------!
585
586
587	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
588	!!$if (itestf.gt.0) then
589	!!$ write(6,*) 'dans flottab apres DETERMIN_TACHE asymetrie sur H pour time=',time
590	!!$ stop
591	!!$else
592	!!$ write(6,*) 'dans flottab apres DETERMIN_TACHE pas d asymetrie sur H pour time=',time
593	!!$
594	!!$end if
595
596	!ice(:,:)=0 Attention, voir si ca marche toujours pour l'Antarctique et heminord !
597
598	!On compte comme englacé uniquement les calottes dont une partie est posée
599
600	do i=3,nx-2
601	do j=3,ny-2
602	test1: if (.not.iceberg(table_out(i,j))) then ! on est pas sur un iceberg
603	if (nb_pts_tache(table_out(i,j)).ge.1) then
604	ice(i,j)=1
605	if (nb_pts_tache(table_out(i,j)).le.10) then ! les petites iles sont en sia
606	! write(6,*) 'petite ile ',i,j
607	flgzmx(i,j)=.false.
608	flgzmx(i+1,j)=.false.
609	flgzmy(i,j)=.false.
610	flgzmy(i,j+1)=.false.
611	gzmx(i:i+1,j)=.false.
612	gzmy(i,j:j+1)=.false.
613	endif
614
615
616	! ici on est sur une tache non iceberg >= 5 points
617	! on teste si la tache n'est pas completement ice stream
618
619	test2: if (icetrim(table_out(i,j))) then ! on a une ile d'ice stream
620
621	mask_tache_ij(:,:)=.false.
622	mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j)) ! pour toute la tache
623
624	smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:))
625	smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:))
626	smax_i=smax_coord(1)
627	smax_j=smax_coord(2)
628
629	!!$ smax_i=0 ; smax_j=0 ; smax_=sealevel
630	!!$ do ii=3,nx-2
631	!!$ do jj=3,ny-2
632	!!$ if (table_out(ii,jj).eq.table_out(i,j)) then
633	!!$ if (s(ii,jj).gt.smax_) then
634	!!$ smax_ =s(ii,jj)
635	!!$ smax_i=ii
636	!!$ smax_j=jj
637	!!$ endif
638	!!$ endif
639	!!$ end do
640	!!$ end do
641
642
643	gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false.
644	gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false.
645	flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false.
646	flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false.
647
648	if (Smax_.le.sealevel) then
649	write(num_tracebug,*)'Attention, une ile avec la surface sous l eau'
650	write(num_tracebug,*)'time=',time,' coord:',smax_i,smax_j
651	end if
652
653	endif test2
654	end if ! endif deplace
655
656	else ! on est sur un iceberg ! test1
657	ice(i,j)=0
658	h(i,j)=1.
659	surnet=H(i,j)*(1.-ro/row)
660	S(i,j)=surnet+sealevel
661	B(i,j)=S(i,j)-H(i,j)
662
663	endif test1
664
665
666	end do
667	end do
668
669
670	!----------------------------------------------
671	! On caracterise le front des ice shelfs/streams
672
673	! call DETERMIN_FRONT
674
675	!----------------------------------------------
676	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
677	!!$if (itestf.gt.0) then
678	!!$ write(6,*) 'dans flottab apres DETERMIN_front asymetrie sur H pour time=',time
679	!!$ stop
680	!!$else
681	!!$ write(6,*) 'dans flottab apres DETERMIN_front pas d asymetrie sur H pour time=',time
682	!!$
683	!!$end if
684
685	endif synchro
686
687	! correction momentanée pour symetrie Heino
688	!where ((.not.flot(:,:)).and.(ice(:,:).eq.0)) H(:,:)=0.
689
690	!fin de routine flottab2
691	!print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j)
692
693
694
695	end subroutine flottab
696	!--------------------------------------------------------------------
697
698	!> SUBROUTINE: determin_tache
699	!! Routine pour la dtermination du numero de tache a effectuer
700	!>
701	subroutine determin_tache
702
703	implicit none
704	integer :: indice
705	integer :: label ! no des taches rencontrées dans le mask
706	integer :: label_max ! no temporaire maxi de tache rencontrées
707	! integer :: mask_nb = 4
708	integer,parameter :: mask_nb = 2 ! version ou on ne compte pas les diagonales
709	integer :: vartemp ! variable temporaire pour reorganiser compt
710	! integer,dimension(mask_nb) :: mask
711	integer,dimension(mask_nb) :: mask
712
713
714	! 1-initialisation
715	!-----------------
716	label_max=1 ! numero de la tache, la premiere tache est notée 1
717	label=1
718	do i=1,n_ta_max
719	compt(i)=i
720	enddo
721	! table_in = .false.
722
723	table_out(:,:) = 0
724	iceberg(:) = .true.
725	icetrim (:) = .true.
726	nb_pts_tache(:) = 0
727
728	! open(unit=100,file="tache.data",status='replace')
729
730	! 2-reperage des taches
731	!----------------------
732	do i=2,nx-1
733	do j=2,ny-1
734
735
736
737	IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------!
738
739	if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then !masque de 2 cases adjacentes
740	! un des voisins est deja en glace
741	mask(1) = table_out(i-1,j)
742	mask(2) = table_out(i,j-1)
743	label = label_max
744
745	!on determine la valeur de la tache minimun (>0) presente ds le masque
746	do indice=1,mask_nb
747	if (mask(indice).gt.0) label=min(label,mask(indice))
748	enddo
749
750	!on fixe la valeur de la tache voisine minimun au point etudie (via label)
751	table_out(i,j)=label
752	!si ce noeud est posé, alors la tache n'est pas un iceberg et iceberg=.F.
753	if (.not.FLOT(I,J)) then
754	iceberg(label)=.false.
755	endif
756
757	!si ce noeud est posé, alors la tache n'est pas un ice stream et icestrim=.F.
758	if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. &
759	(.not.gzmy(i,j).and..not.gzmy(i,j+1))) then
760	icetrim(label)=.false.
761	endif
762
763	! si 2 taches differentes sont dans le masque, il faut les identifier dans compt
764	! i.e. les plus grands numeros correspondent au plus petit
765	! on lui affecte le numero de la tache fondamentale avec un signe -
766	! pour indiquer le changement
767
768	do indice=1,mask_nb
769	if(mask(indice).gt.label) then
770	compt(mask(indice))=-label
771	endif
772	enddo
773
774	else !aucun des voisins est une tache
775	table_out(i,j)= label_max
776	compt(label_max)=label_max
777	!si ce noeud est posé, alors la ache n'est pas un iceberg et iceberg=.F.
778	if (.not.FLOT(I,J)) then
779	iceberg(label_max)=.false.
780	endif
781
782	!si ce noeud est posé, alors le tache n'est pas un ice stream et icestrim=.F.
783	if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. &
784	(.not.gzmy(i,j).and..not.gzmy(i,j+1))) then
785	icetrim(label)=.false.
786	endif
787
788	label_max = label_max+1
789	if (label_max.gt.n_ta_max) print*,'trop de taches=',label_max
790	endif
791
792
793	else !on est pas sur une tache----------------------------------------------
794	table_out(i,j)=0 ! Pas necessaire (reecrit 0 sur 0)
795	endif !---------------------------------------------------------------------
796
797
798	enddo
799	enddo
800
801
802
803	! On reorganise compt en ecrivant le numero de la tache fondamentale
804	! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité)
805	! On indique aussi le nb de point que contient chaque taches (nb_pts_tache)
806	do indice=1,label_max
807	vartemp = compt(indice)
808	if (compt(indice).lt.0) then
809	compt(indice)= compt(-vartemp)
810	if (.not.iceberg(indice)) iceberg(-vartemp)=.false.
811	if (.not.icetrim(indice)) icetrim(-vartemp)=.false.
812	endif
813	enddo
814
815	do i=1,nx
816	do j=1,ny
817	if (table_out(i,j).ne.0) then
818	table_out(i,j)=compt(table_out(i,j))
819	nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1
820	endif
821	enddo
822	enddo
823
824
825
826
827	!!$tablebis(:,:)=table_out(:,:)
828	!!$do j=1,ny
829	!!$ do i=1,nx
830	!!$ if (tablebis(i,j).ne.0) then ! tache de glace
831	!!$ numtache=table_out(i,j)
832	!!$ nb_pt=count(table_out(:,:).eq.numtache) ! compte tous les points de la tache
833	!!$ nb_pts_tache(table_out(i,j))=nb_pt !
834	!!$
835	!!$ where (table_out(:,:).eq.numtache)
836	!!$ tablebis(:,:)=0 ! la table de tache est remise a 0 pour eviter de repasser
837	!!$ end where
838	!!$ write(6,*) i,j,nb_pt,table_out(i,j)
839	!!$ endif
840	!!$ enddo
841	!!$enddo
842
843	!!$do j=1,ny
844	!!$ do i=1,nx
845	!!$ debug_3d(i,j,56)=nb_pts_tache(table_out(i,j))
846	!!$ end do
847	!!$end do
848	!!$
849	end subroutine determin_tache
850	!----------------------------------------------------------------------
851	!> SUBROUTINE: determin_front
852	!!Routine pour la determination du front
853	!>
854	subroutine determin_front
855
856	integer :: i_moins1,i_plus1,i_plus2
857	integer :: j_moins1,j_plus1,j_plus2
858
859	do i=3,nx-2
860	do j=3,ny-2
861
862	surice:if (ice(i,j).eq.0) then
863	do pmx=-1,1,2
864	do pmy=-1,1,2
865
866	diagice : if (ice(i+pmx,j+pmy).eq.1) then
867
868	if ((ice(i+pmx,j)+ice(i,j+pmy).eq.2)) then ! test (i) pour eviter les langues
869	! de glaces diagonales en coin(26dec04)
870	if ((ice(i+2pmx,j).eq.1.and.ice(i+2pmx,j+pmy).eq.0).or.&
871	(ice(i,j+2pmy).eq.1.and.ice(i+pmx,j+2pmy).eq.0)) then
872	ice(i,j)=1
873	h(i,j)=max(1.,h(i,j))
874	endif
875
876	! test (i) pour eviter les langues de glaces diagonales :
877	! mouvement du cheval aux echecs
878
879	if ((ice(i+2pmx,j+pmy)+ice(i+pmx,j+2pmy).eq.1)) then
880	if (ice(i+2*pmx,j+pmy).eq.1.and. &
881	(ice(i+2pmx,j+2pmy)+ice(i,j+2*pmy)).ge.1) then
882	ice(i,j)=1
883	h(i,j)=max(1.,h(i,j))
884	endif
885	if (ice(i+pmx,j+2*pmy).eq.1.and. &
886	(ice(i+2pmx,j+2pmy)+ice(i+2*pmx,j)).ge.1) then
887	ice(i,j)=1
888	h(i,j)=max(1.,h(i,j))
889	endif
890
891	! test (ii) pour eviter les langues de glaces diagonales :
892	! le point glace ice(i+pmx,j+pmy) a :
893	! - ses 4 voisins frontaux en glace
894	! - mais 2 voisins vides diagonalement opposes
895
896	elseif ((ice(i+2pmx,j+pmy)+ice(i+pmx,j+2pmy).eq.2) &
897	.and.ice(i+2pmx,j+2pmy).eq.0) then
898
899	! test (iii) pour faire les tests (i) et (ii)
900	ice(i,j)=1
901	h(i,j)=max(1.,h(i,j))
902	ice(i+2pmx,j+2pmy)=1
903	h(i+2pmx,j+2pmy)=max(1.,h(i+2pmx,j+2pmy))
904	endif
905	endif
906	endif diagice
907	enddo
908	enddo
909	endif surice
910	end do
911	end do
912
913
914	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
915	!!$if (itestf.gt.0) then
916	!!$ write(6,*) 'dans front avant remplissage baies asymetrie sur H pour time=',time
917	!!$ stop
918	!!$else
919	!!$ write(6,*) 'dans front avant remplissage baies pas d asymetrie sur H pour time=',time
920	!!$
921	!!$end if
922
923
924	! print*,'dans remplissage baies',time
925	baies: do k=1,2
926	do j=1,ny
927	do i=1,nx
928
929	surice_xy: if (ice(i,j).eq.0) then
930	i_moins1=max(i-1,1)
931	j_moins1=max(j-1,1)
932	i_plus1=min(i+1,nx)
933	j_plus1=min(j+1,ny)
934	i_plus2=min(i+2,nx)
935	j_plus2=min(j+2,ny)
936
937	! test (iii) pour trouver les baies de largeur 1 ou 2 cases
938	! et combler les trous si ce sont des baies
939	! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles
940	! baies horizontales
941
942	if (ice(i_moins1,j).eq.1.and.(ice(i_plus1,j).eq.1.or.ice(i_plus2,j).eq.1)) then
943	if (ice(i,j_moins1).eq.1.or.ice(i,j_plus1).eq.1) then ! ice(i,j)=1
944	ice(i,j)=1
945	H(i,j)=max(1.,H(i,j))
946	endif
947	endif
948
949
950	if (ice(i,j_moins1).eq.1.and.(ice(i,j_plus1).eq.1.or.ice(i,j_plus2).eq.1)) then
951	if (ice(i_moins1,j).eq.1.or.ice(i_plus1,j).eq.1) then !ice(i,j)=1
952	ice(i,j)=1
953	H(i,j)=max(1.,H(i,j))
954	endif
955	endif
956
957	endif surice_xy
958	end do
959	end do
960	end do baies
961
962	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
963	!!$if (itestf.gt.0) then
964	!!$ write(6,*) 'dans front apres remplissage baies asymetrie sur H pour time=',time
965	!!$ stop
966	!!$else
967	!!$ write(6,*) 'dans front apres remplissage baies pas d asymetrie sur H pour time=',time
968	!!$
969	!!$end if
970
971
972	do i=2,nx-1
973	do j=2,ny-1
974
975	if (ice(i,j).eq.1) then ! test si ice=1
976
977	! if ice, on determine front...
978	! ainsi, front=0 sur les zones = 0
979
980	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
981	!front= le nb de faces en contact avec un voisin englacé
982	endif
983	end do
984	end do
985
986	! traitement des bords. On considere que l'exterieur n'a pas de glace
987	! attention ce n'est vrai que sur la grande grille
988
989
990	do j=2,ny-1
991	i=1
992	front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
993	i=nx
994	front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
995	end do
996
997	do i=2,nx-1
998	j=1
999	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))
1000	j=ny
1001	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
1002	end do
1003
1004	! traitement des coins
1005
1006	front(1,1)=ice(2,1)+ice(2,1)
1007	front(1,ny)=ice(2,ny)+ice(1,ny-1)
1008	front(nx,1)=ice(nx,2)+ice(nx-1,1)
1009	front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)
1010
1011	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
1012	!!$if (itestf.gt.0) then
1013	!!$ write(6,*) 'dans front apres front asymetrie sur H pour time=',time
1014	!!$ stop
1015	!!$else
1016	!!$ write(6,*) 'dans front apres front pas d asymetrie sur H pour time=',time
1017	!!$
1018	!!$end if
1019
1020	! on ne compte pas les taches de glace de 2 cases (horizontales ou verticales)
1021	! en fait, si ces deux cases sont flottantes, il faut enlever les icebergs
1022	! de n'importe quelle taille).
1023	! si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites
1024
1025
1026	do j=1,ny
1027	do i=1,nx-1
1028	if (front(i,j).eq.1) then
1029	if (front(i+1,j).eq.1) then
1030	ice(i,j)=0
1031	ice(i+1,j)=0
1032	front(i,j)=0
1033	front(i+1,j)=0
1034	endif
1035	endif
1036	end do
1037	end do
1038
1039	do j=1,ny-1
1040	do i=1,nx
1041	if (front(i,j).eq.1) then
1042	if (front(i,j+1).eq.1) then
1043	ice(i,j)=0
1044	ice(i,j+1)=0
1045	front(i,j)=0
1046	front(i,j+1)=0
1047	endif
1048	end if
1049	end do
1050	end do
1051
1052	!isolx signifie pas de voisins en x
1053	!isoly signifie pas de voisins en y
1054	!remarque :
1055	!si isolx/y=.true. alors frontfacex/y=0 (a la fois +1 & -1 or +1-1=0)
1056
1057	! calcul de frontfacex et isolx
1058	do j=1,ny
1059	do i=2,nx-1
1060
1061	if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3)
1062
1063	if ((ice(i-1,j)+ice(i+1,j)).lt.2) then ! il y a un front // a x
1064
1065	if ((ice(i-1,j)+ice(i+1,j)).eq.0) then
1066	isolx(i,j)=.true.
1067	elseif (ice(i-1,j).eq.0) then
1068	frontfacex(i,j)=-1 ! front i-1 \|i i+1
1069	else
1070	frontfacex(i,j)=+1 ! front i-1 i\| i+1
1071	endif
1072	endif
1073	end if !fin du test il y a un front
1074
1075	end do
1076	end do
1077
1078	! calcul de frontfacey et isoly
1079	do j=2,ny-1
1080	do i=1,nx
1081
1082	if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3)
1083
1084	if ((ice(i,j-1)+ice(i,j+1)).lt.2) then ! il y a un front // a y
1085
1086	if ((ice(i,j-1)+ice(i,j+1)).eq.0) then
1087	isoly(i,j)=.true. !front j-1 \|j\| j+1
1088	elseif (ice(i,j-1).eq.0) then
1089	frontfacey(i,j)=-1 !front j-1 \|j j+1
1090	else
1091	frontfacey(i,j)=+1 !front j-1 j\| j+1
1092	endif
1093	endif
1094	end if !fin du test il y a un front
1095
1096	end do
1097	end do
1098
1099
1100
1101	! traitement des bords. On considere que l'exterieur n'a pas de glace
1102	! attention ce n'est vrai que sur la grande grille
1103
1104
1105	do j=2,ny-1
1106	i=1
1107	if (front(i,j).ge.1) then
1108	if (ice(i+1,j).eq.0) then
1109	isolx(i,j)=.true.
1110	else
1111	frontfacex(i,j)=-1
1112	endif
1113	end if
1114	i=nx
1115	if (front(i,j).ge.1) then
1116	if (ice(i-1,j).eq.0) then
1117	isolx(i,j)=.true.
1118	else
1119	frontfacex(i,j)=1
1120	endif
1121	end if
1122	end do
1123
1124	do i=2,nx-1
1125	j=1
1126	if (front(i,j).ge.1) then
1127	if (ice(i,j+1).eq.0) then
1128	isoly(i,j)=.true.
1129	else
1130	frontfacey(i,j)=-1
1131	endif
1132	end if
1133	j=ny
1134	if (front(i,j).ge.1) then
1135	if (ice(i,j-1).eq.0) then
1136	isoly(i,j)=.true.
1137	else
1138	frontfacey(i,j)=1
1139	endif
1140	end if
1141	end do
1142
1143
1144
1145	return
1146	end subroutine determin_front
1147	!------------------------------------------------------------------------------
1148
1149	end module flottab_mod
1150
1151

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

Download in other formats: