Context Navigation

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

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

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

Download in other formats: