Context Navigation

source: branches/GRISLIv3/SOURCES/flottab2-0.7.f90 @ 468

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

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

Download in other formats: