Context Navigation

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

Last change on this file since 71 was 71, checked in by dumas, 8 years ago
First parallelization instructions with OpenMP tested in debug : exactly the same results \| isostasy called only every 50 years
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)
633	!$OMP DO
634	do i=3,nx-2
635	do j=3,ny-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	!!$OMP PARALLEL PRIVATE(mask,label,indice)
784	!!$OMP DO
785	do i=2,nx-1
786	do j=2,ny-1
787
788
789
790	IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------!
791
792	if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then !masque de 2 cases adjacentes
793	! un des voisins est deja en glace
794	mask(1) = table_out(i-1,j)
795	mask(2) = table_out(i,j-1)
796	label = label_max
797
798	!on determine la valeur de la tache minimun (>0) presente ds le masque
799	do indice=1,mask_nb
800	if (mask(indice).gt.0) label=min(label,mask(indice))
801	enddo
802	!cdc label=min(label,minval(mask(:), mask=mask > 0))
803
804	!on fixe la valeur de la tache voisine minimun au point etudie (via label)
805	table_out(i,j)=label
806	!si ce noeud est posé, alors la tache n'est pas un iceberg et iceberg=.F.
807	if (.not.FLOT(I,J)) then
808	iceberg(label)=.false.
809	endif
810
811	!si ce noeud est posé, alors la tache n'est pas un ice stream et icestrim=.F.
812	if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. &
813	(.not.gzmy(i,j).and..not.gzmy(i,j+1))) then
814	icetrim(label)=.false.
815	endif
816
817	! si 2 taches differentes sont dans le masque, il faut les identifier dans compt
818	! i.e. les plus grands numeros correspondent au plus petit
819	! on lui affecte le numero de la tache fondamentale avec un signe -
820	! pour indiquer le changement
821
822	do indice=1,mask_nb
823	if(mask(indice).gt.label) then
824	compt(mask(indice))=-label
825	endif
826	enddo
827
828	else !aucun des voisins est une tache
829	table_out(i,j)= label_max
830	compt(label_max)=label_max
831	!si ce noeud est posé, alors la ache n'est pas un iceberg et iceberg=.F.
832	if (.not.FLOT(I,J)) then
833	iceberg(label_max)=.false.
834	endif
835
836	!si ce noeud est posé, alors le tache n'est pas un ice stream et icestrim=.F.
837	if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and. &
838	(.not.gzmy(i,j).and..not.gzmy(i,j+1))) then
839	icetrim(label)=.false.
840	endif
841
842	label_max = label_max+1
843	if (label_max.gt.n_ta_max) print*,'trop de taches=',label_max
844	endif
845
846
847	else !on est pas sur une tache----------------------------------------------
848	table_out(i,j)=0 ! Pas necessaire (reecrit 0 sur 0)
849	endif !---------------------------------------------------------------------
850
851
852	enddo
853	enddo
854	!!$OMP END DO
855	!!$OMP END PARALLEL
856
857
858	! On reorganise compt en ecrivant le numero de la tache fondamentale
859	! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité)
860	! On indique aussi le nb de point que contient chaque taches (nb_pts_tache)
861
862	do indice=1,label_max
863	vartemp = compt(indice)
864	if (compt(indice).lt.0) then
865	compt(indice)= compt(-vartemp)
866	if (.not.iceberg(indice)) iceberg(-vartemp)=.false.
867	if (.not.icetrim(indice)) icetrim(-vartemp)=.false.
868	endif
869	enddo
870
871	!!$OMP PARALLEL
872	!!$OMP DO REDUCTION(+:nb_pts_tache)
873	do i=1,nx
874	do j=1,ny
875	if (table_out(i,j).ne.0) then
876	table_out(i,j)=compt(table_out(i,j))
877	nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1
878	endif
879	enddo
880	enddo
881	!!$OMP END DO
882	!!$OMP END PARALLEL
883
884
885	!!$tablebis(:,:)=table_out(:,:)
886	!!$do j=1,ny
887	!!$ do i=1,nx
888	!!$ if (tablebis(i,j).ne.0) then ! tache de glace
889	!!$ numtache=table_out(i,j)
890	!!$ nb_pt=count(table_out(:,:).eq.numtache) ! compte tous les points de la tache
891	!!$ nb_pts_tache(table_out(i,j))=nb_pt !
892	!!$
893	!!$ where (table_out(:,:).eq.numtache)
894	!!$ tablebis(:,:)=0 ! la table de tache est remise a 0 pour eviter de repasser
895	!!$ end where
896	!!$ write(6,*) i,j,nb_pt,table_out(i,j)
897	!!$ endif
898	!!$ enddo
899	!!$enddo
900
901	!!$do j=1,ny
902	!!$ do i=1,nx
903	!!$ debug_3d(i,j,56)=nb_pts_tache(table_out(i,j))
904	!!$ end do
905	!!$end do
906	!!$
907	end subroutine determin_tache
908	!----------------------------------------------------------------------
909	!> SUBROUTINE: determin_front
910	!!Routine pour la determination du front
911	!>
912	subroutine determin_front
913	!!$ USE OMP_LIB
914	integer :: i_moins1,i_plus1,i_plus2
915	integer :: j_moins1,j_plus1,j_plus2
916
917	!!$OMP PARALLEL
918	!!$OMP DO
919	do j=3,ny-2
920	do i=3,nx-2
921
922	surice:if (ice(i,j).eq.0) then
923	do pmx=-1,1,2
924	do pmy=-1,1,2
925
926	diagice : if (ice(i+pmx,j+pmy).eq.1) then
927
928	if ((ice(i+pmx,j)+ice(i,j+pmy).eq.2)) then ! test (i) pour eviter les langues
929	! de glaces diagonales en coin(26dec04)
930	if ((ice(i+2pmx,j).eq.1.and.ice(i+2pmx,j+pmy).eq.0).or.&
931	(ice(i,j+2pmy).eq.1.and.ice(i+pmx,j+2pmy).eq.0)) then
932	ice(i,j)=1
933	h(i,j)=max(1.,h(i,j))
934	endif
935
936	! test (i) pour eviter les langues de glaces diagonales :
937	! mouvement du cheval aux echecs
938
939	if ((ice(i+2pmx,j+pmy)+ice(i+pmx,j+2pmy).eq.1)) then
940	if (ice(i+2*pmx,j+pmy).eq.1.and. &
941	(ice(i+2pmx,j+2pmy)+ice(i,j+2*pmy)).ge.1) then
942	ice(i,j)=1
943	h(i,j)=max(1.,h(i,j))
944	endif
945	if (ice(i+pmx,j+2*pmy).eq.1.and. &
946	(ice(i+2pmx,j+2pmy)+ice(i+2*pmx,j)).ge.1) then
947	ice(i,j)=1
948	h(i,j)=max(1.,h(i,j))
949	endif
950
951	! test (ii) pour eviter les langues de glaces diagonales :
952	! le point glace ice(i+pmx,j+pmy) a :
953	! - ses 4 voisins frontaux en glace
954	! - mais 2 voisins vides diagonalement opposes
955
956	elseif ((ice(i+2pmx,j+pmy)+ice(i+pmx,j+2pmy).eq.2) &
957	.and.ice(i+2pmx,j+2pmy).eq.0) then
958
959	! test (iii) pour faire les tests (i) et (ii)
960	ice(i,j)=1
961	h(i,j)=max(1.,h(i,j))
962	ice(i+2pmx,j+2pmy)=1
963	h(i+2pmx,j+2pmy)=max(1.,h(i+2pmx,j+2pmy))
964	endif
965	endif
966	endif diagice
967	enddo
968	enddo
969	endif surice
970	end do
971	end do
972	!!$OMP END DO
973
974	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
975	!!$if (itestf.gt.0) then
976	!!$ write(6,*) 'dans front avant remplissage baies asymetrie sur H pour time=',time
977	!!$ stop
978	!!$else
979	!!$ write(6,*) 'dans front avant remplissage baies pas d asymetrie sur H pour time=',time
980	!!$
981	!!$end if
982
983
984	! print*,'dans remplissage baies',time
985
986	baies: do k=1,2
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	end do baies
1024
1025	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
1026	!!$if (itestf.gt.0) then
1027	!!$ write(6,*) 'dans front apres remplissage baies asymetrie sur H pour time=',time
1028	!!$ stop
1029	!!$else
1030	!!$ write(6,*) 'dans front apres remplissage baies pas d asymetrie sur H pour time=',time
1031	!!$
1032	!!$end if
1033
1034	!!$OMP DO
1035	do j=2,ny-1
1036	do i=2,nx-1
1037
1038	if (ice(i,j).eq.1) then ! test si ice=1
1039
1040	! if ice, on determine front...
1041	! ainsi, front=0 sur les zones = 0
1042
1043	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
1044	!front= le nb de faces en contact avec un voisin englacé
1045	endif
1046	end do
1047	end do
1048	!!$OMP END DO
1049
1050	! traitement des bords. On considere que l'exterieur n'a pas de glace
1051	! attention ce n'est vrai que sur la grande grille
1052
1053	!!$OMP DO PRIVATE(i)
1054	do j=2,ny-1
1055	i=1
1056	front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
1057	i=nx
1058	front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
1059	end do
1060	!!$OMP END DO
1061
1062	!!$OMP DO PRIVATE(j)
1063	do i=2,nx-1
1064	j=1
1065	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))
1066	j=ny
1067	front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
1068	end do
1069	!!$OMP END DO
1070	! traitement des coins
1071
1072	front(1,1)=ice(2,1)+ice(2,1)
1073	front(1,ny)=ice(2,ny)+ice(1,ny-1)
1074	front(nx,1)=ice(nx,2)+ice(nx-1,1)
1075	front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)
1076
1077	!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
1078	!!$if (itestf.gt.0) then
1079	!!$ write(6,*) 'dans front apres front asymetrie sur H pour time=',time
1080	!!$ stop
1081	!!$else
1082	!!$ write(6,*) 'dans front apres front pas d asymetrie sur H pour time=',time
1083	!!$
1084	!!$end if
1085
1086	! on ne compte pas les taches de glace de 2 cases (horizontales ou verticales)
1087	! en fait, si ces deux cases sont flottantes, il faut enlever les icebergs
1088	! de n'importe quelle taille).
1089	! si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites
1090
1091	!!$OMP DO
1092	do j=1,ny
1093	do i=1,nx-1
1094	if (front(i,j).eq.1) then
1095	if (front(i+1,j).eq.1) then
1096	ice(i,j)=0
1097	ice(i+1,j)=0
1098	front(i,j)=0
1099	front(i+1,j)=0
1100	endif
1101	endif
1102	end do
1103	end do
1104	!!$OMP END DO
1105
1106	!!$OMP DO
1107	do j=1,ny-1
1108	do i=1,nx
1109	if (front(i,j).eq.1) then
1110	if (front(i,j+1).eq.1) then
1111	ice(i,j)=0
1112	ice(i,j+1)=0
1113	front(i,j)=0
1114	front(i,j+1)=0
1115	endif
1116	end if
1117	end do
1118	end do
1119	!!$OMP END DO
1120
1121	!isolx signifie pas de voisins en x
1122	!isoly signifie pas de voisins en y
1123	!remarque :
1124	!si isolx/y=.true. alors frontfacex/y=0 (a la fois +1 & -1 or +1-1=0)
1125
1126	! calcul de frontfacex et isolx
1127	!!$OMP DO
1128	do j=1,ny
1129	do i=2,nx-1
1130
1131	if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3)
1132
1133	if ((ice(i-1,j)+ice(i+1,j)).lt.2) then ! il y a un front // a x
1134
1135	if ((ice(i-1,j)+ice(i+1,j)).eq.0) then
1136	isolx(i,j)=.true.
1137	elseif (ice(i-1,j).eq.0) then
1138	frontfacex(i,j)=-1 ! front i-1 \|i i+1
1139	else
1140	frontfacex(i,j)=+1 ! front i-1 i\| i+1
1141	endif
1142	endif
1143	end if !fin du test il y a un front
1144
1145	end do
1146	end do
1147	!!$OMP END DO
1148
1149	! calcul de frontfacey et isoly
1150	!!$OMP DO
1151	do j=2,ny-1
1152	do i=1,nx
1153
1154	if (front(i,j).ge.1.and.front(i,j).le.3) then !front(entre 1 et 3)
1155
1156	if ((ice(i,j-1)+ice(i,j+1)).lt.2) then ! il y a un front // a y
1157
1158	if ((ice(i,j-1)+ice(i,j+1)).eq.0) then
1159	isoly(i,j)=.true. !front j-1 \|j\| j+1
1160	elseif (ice(i,j-1).eq.0) then
1161	frontfacey(i,j)=-1 !front j-1 \|j j+1
1162	else
1163	frontfacey(i,j)=+1 !front j-1 j\| j+1
1164	endif
1165	endif
1166	end if !fin du test il y a un front
1167
1168	end do
1169	end do
1170	!!$OMP END DO
1171
1172
1173	! traitement des bords. On considere que l'exterieur n'a pas de glace
1174	! attention ce n'est vrai que sur la grande grille
1175
1176	!!$OMP DO PRIVATE(i)
1177	do j=2,ny-1
1178	i=1
1179	if (front(i,j).ge.1) then
1180	if (ice(i+1,j).eq.0) then
1181	isolx(i,j)=.true.
1182	else
1183	frontfacex(i,j)=-1
1184	endif
1185	end if
1186	i=nx
1187	if (front(i,j).ge.1) then
1188	if (ice(i-1,j).eq.0) then
1189	isolx(i,j)=.true.
1190	else
1191	frontfacex(i,j)=1
1192	endif
1193	end if
1194	end do
1195	!!$OMP END DO
1196
1197	!!$OMP DO PRIVATE(j)
1198	do i=2,nx-1
1199	j=1
1200	if (front(i,j).ge.1) then
1201	if (ice(i,j+1).eq.0) then
1202	isoly(i,j)=.true.
1203	else
1204	frontfacey(i,j)=-1
1205	endif
1206	end if
1207	j=ny
1208	if (front(i,j).ge.1) then
1209	if (ice(i,j-1).eq.0) then
1210	isoly(i,j)=.true.
1211	else
1212	frontfacey(i,j)=1
1213	endif
1214	end if
1215	end do
1216	!!$OMP END DO
1217	!!OMP END PARALLEL
1218
1219	return
1220	end subroutine determin_front
1221	!------------------------------------------------------------------------------
1222
1223	end module flottab_mod
1224
1225

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

Download in other formats: