Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

remplimat-ant-0.5-40km.f90 @ 159

Last change on this file since 159 was 62, checked in by dumas, 8 years ago
Move unused sources files in Old-sources directory
File size: 75.9 KB

Line
1	!> \file remplimat-ant-0.5-40km.f90
2	!! TOOOOOOO DOOOOOO
3	!<
4	!====================================================================!
5	SUBROUTINE REMPLIDOM(N1,N2,uxprec,uyprec,uxnew,uynew,ifail) !
6	! !
7	! Vince 17 Jan !
8	! Calcul de la vitesse !
9	! des coins -> 25 fev !
10	! !
11	! Modification du schema numerique 3 Mars 95 (stabilite) !
12	! Aout 95 (optimisation) !
13	! Modifications Catherine Mars 97 !
14	! gradient conjugue utilisant blas !
15	! !
16	! Modification catherine Avril 97 !
17	! DECOMPOSITION DE DOMAINE !
18	! Catherine Mai 97 application a l\'Antarctique !
19	! !
20	! Reecriture Vincent novembre 2003 !
21	! On remanie les conditions aux limites !
22	! ------------------------------------------------------ !
23	! numdom numero du domaine traite !
24	! id1,id2 bornes selon x !
25	! jd1,jd2 bornes selon y !
26	! !
27	! uxprex(n1,n2) !
28	! uyprec(n1,n2) vitesses de l\'iteration precedente !
29	! !
30	! uxnew(n1,n2) !
31	! uynew(n1,n2) uynew resultat de cette iteration !
32	! !
33	!====================================================================!
34
35	USE module3D_phy
36	USE EQ_ELLIPTIQUE_MOD
37
38	IMPLICIT NONE
39
40	! dimensionnement momentane de slv
41	real, dimension(nx,ny) :: slv
42
43	! 1_Declaration des variables
44	! ----------------------------
45
46	integer nblig
47	integer aaai,aaaj
48	real aaa1
49	integer iinf,isup,jinf,jsup
50	integer iii,jjj ! position du noeud de l'echec de sgbsv
51	real aaa
52	INTEGER LDB,NRHS,LDAB,N1,N2
53	INTEGER kzero,knonul,ITER
54	LOGICAL CONJGR!,IREORD
55	INTEGER nb_err ! nombre de passage dans sgbsv
56	! apres le 3eme passage infructueux on stoppoe
57	!!!!!! INTEGER IPIV(NPTMAX)
58
59
60	INTEGER CTEURLIG,SNR,INEARX
61	INTEGER IFLOTMX(NX,NY),IFLOTMY(NX,NY),ifail
62
63	! dimensions liees au domaine
64	! INTEGER NUMDOM,ID1,ID2,JD1,JD2
65	REAL UXPREC(NX,NY),UYPREC(NX,NY),UXNEW(NX,NY),UYNEW(NX,NY)
66	!
67	!Pour le remplissage de la matice
68	REAL SCAL
69	INTEGER IIMAT
70	!INTEGER,DIMENSION(NX,NY) :: INEARX ! la simple variable est devenu
71	! tableau !!!!!!!! voir si on peut la remplacer par un tableau deja
72	!!!!!!!!!existant car elle est utile que temporairement
73	!!!!!!!!! dans §2_Initialisation des variables
74	REAL,DIMENSION(NX,NY) :: TUMIMJ, TUMIJ, TUMIPJ, &
75	TUIMJ , TUIJ , TUIPJ , &
76	TUPIMJ, TUPIJ, TUPIPJ, &
77	TVMIJ,TVMIPJ,TVIJ,TVIPJ !Pour le calcul des U
78
79	REAL,DIMENSION(NX,NY) :: SVMIMJ, SVMIJ, SVMIPJ, &
80	SVIMJ , SVIJ , SVIPJ , &
81	SVPIMJ, SVPIJ, SVPIPJ, &
82	SUIMJ,SUIJ,SUPIMJ,SUPIJ !Pour le calcul des V
83
84	REAL,DIMENSION(NX,NY) :: OPPOSX,OPPOSY!Pression laterale contrebalancant le mouvement pour u et v
85	logical,DIMENSION(NX,NY) ::okumatest,okvmatest ! test de controle pour
86
87	! ligu, ligv transporte depuis 3D
88	integer, dimension(NX,NY) :: LIGU ! numero de ligne de U dans remplidom
89	integer, dimension(NX,NY) :: LIGV ! numero de ligne de V dans remplidom
90
91	real upr,unw, emmax
92	real moteur,beta!,unorm
93	INTEGER,DIMENSION(:,:),allocatable:: posligu, posligv !transpose de ligu,ligv
94	INTEGER,DIMENSION(:),allocatable:: ipiv !7eme parametre de sgbsv (pour permutations)
95	CHARACTER(len=8),DIMENSION(NX,NY,2) :: METHOD
96	INTEGER,DIMENSION(NX,NY,2) :: METHOD2
97
98	integer :: num_err = 266
99
100
101
102	! sorties debug
103	open(num_err,file='erreurs-remplimat')
104	write(num_err,*)' time=',time
105
106	uxnew(:,:)=0.
107	uynew(:,:)=0.
108
109	! attention, pour l'instant slv est local et juste attribue ici
110	slv(:,:)=sealevel
111
112	kzero =0
113	nb_err=0
114
115	! 2_Initialisation des variables
116	! ----------------------------------
117
118	METHOD(:,:,:)='000000'
119	METHOD2(:,:,:)=-2
120	NRHS=1
121	LDAB=2*KL+KU+1
122	LDB=NPTMAX
123	INEARX=1
124	CTEURLIG = 0
125
126	iii=0; jjj=0
127	LIGU(:,:)=0
128	LIGV(:,:)=0
129
130	OKUMAT(:,:)=.FALSE.
131	OKVMAT(:,:)=.FALSE.
132	OKUMATEST(:,:)=.FALSE.! TEST
133	OKVMATEST(:,:)=.FALSE.
134
135	! attribution de iflotmx et iflotmy pour remplacer les operations
136	! sur les logiques, iflot est vrai quand soit flotmx ou gzmx est vrai (soit flgzmx)
137
138
139	do J=1,NY
140	do I=1,NX
141
142
143	if (FLGZMX(I,J)) then
144	IFLOTMX(I,J)=1
145	FROTMX(I,J)=min(abs(betamx(i,j)),betamax)
146	else
147	IFLOTMX(I,J)=0.
148	FROTMX(I,J)=abs(betamx(i,j))
149	endif
150
151
152	if (FLGZMY(I,J)) then
153	IFLOTMY(I,J)=1
154	FROTMY(I,J)= min(abs(betamy(i,j)),betamax)
155	else
156	IFLOTMY(I,J)=0
157	FROTMY(I,J)=abs(betamy(i,j))
158	endif
159
160	end do
161	end do
162
163	!testvincent
164	! print*,'emplimat',time,sealevel,tafor
165	!si le point ou un des voisins et stream/shelf inearx.ge.1
166	!!!!!DO I=ID1,ID2
167	!!!!!!!!DO J=JD1,JD2
168	DO j=2,ny !J=1,NY
169	DO i=2,nx! I=1,NX
170	!!! Faire des tests sur les shifts pour amelioration code
171	!!! plutot somme de sous-tableaux
172	IF (ICE(I,J).EQ.1) THEN
173	INEARX=IFLOTMX(I-1,J) +IFLOTMX(I,J) +IFLOTMX(I+1,J) &
174	+IFLOTMX(I-1,J-1)+IFLOTMX(I,J-1)+IFLOTMX(I+1,J-1) &
175	+IFLOTMX(I-1,J+1)+IFLOTMX(I,J+1)+IFLOTMX(I+1,J+1) &
176	+IFLOTMY(I-1,J) +IFLOTMY(I,J) +IFLOTMY(I+1,J) &
177	+IFLOTMY(I-1,J-1)+IFLOTMY(I,J-1)+IFLOTMY(I+1,J-1) &
178	+IFLOTMY(I-1,J+1)+IFLOTMY(I,J+1)+IFLOTMY(I+1,J+1)
179	ELSE
180	INEARX=0
181	END IF
182
183	!print*,('i=',i,'j=',j,'INEARX=',INEARX)
184	IF (INEARX.ge.1) THEN
185	OKUMAT(I,J) = .TRUE.
186	OKUMAT(I+1,J) = .TRUE. ! redondant ???
187	OKVMAT(I,J) = .TRUE.
188	OKVMAT(I,J+1) = .TRUE.
189	! attention a la methode de calcul, les points Ui+1j sont ejectes aux
190	! calandes
191	ENDIF
192	END DO
193	END DO
194
195	DO I=1,NX
196	DO J=1,NY
197	IF (OKUMAT(I,J)) THEN
198	cteurlig=cteurlig+1
199	LIGU(i,j)=cteurlig
200	ENDIF
201	IF (OKVMAT(I,J)) THEN
202	cteurlig=cteurlig+1
203	LIGV(i,j)=cteurlig
204	ENDIF
205	END DO
206	END DO
207
208	! nblig est le nombre d'equation lineaire a resoudre dans mmat
209	NBLIG = CTEURLIG
210	!!!!!!Coeff pour mmat et sgbsv !!!!!!
211	! print*,'avant largeur de bande, ku=',ku,'kl=',kl
212	call larg_bande(ku,kl)
213	! print*,'apres largeur de bande, ku=',ku,'kl=',kl
214	if (ku.lt.0.or.kl.lt.0) then
215	print*,'ku=',ku,'kl=',kl
216	print*,'pas de shelves a traiter : sortie de remplimat'
217	ifail=0
218	return
219	endif
220	kdc2=ku+kl+1
221	! print*,'kdc2=',kdc2
222	ldab=2*kl+ku+1
223	! print*,'ldab=',ldab
224	ldb=nblig
225
226	!test sur l'ancienne taille de mmat
227	! si mmat est trop petite, on l'agrandie a la taille nessessaire
228	if (nblig.gt.nptmax.or.ldab.gt.LDBMAX) then
229	!print*,'redefinition des matrices',nblig,nptmax,ldab,LDBMAX
230	nptmax=max(nptmax,nblig)
231	LDBMAX=max(LDBMAX,ldab)
232	call DEFINITION_MATRICE(.FALSE.,nptmax,ldbmax)
233	endif
234
235	allocate(POSLIGU(nblig,2),stat=err)
236	if (err/=0) then
237	print *,"Erreur à l'allocation de POSLIGU",err
238	stop 4
239	end if
240	allocate(POSLIGV(nblig,2),stat=err)
241	if (err/=0) then
242	print *,"Erreur à l'allocation de POSLIGV",err
243	stop 4
244	end if
245	allocate(IPIV(nblig),stat=err)
246	if (err/=0) then
247	print *,"Erreur à l'allocation de IPIV",err
248	stop 4
249	end if
250
251	cteurlig=0
252	POSLIGU(:,:)=0
253	POSLIGV(:,:)=0
254
255	DO I=1,NX
256	DO J=1,NY
257	IF (OKUMAT(I,J)) THEN
258	cteurlig=cteurlig+1
259	LIGU(i,j)=cteurlig
260	POSLIGU(cteurlig,1)=i
261	POSLIGU(cteurlig,2)=j
262	ENDIF
263	IF (OKVMAT(I,J)) THEN
264	cteurlig=cteurlig+1
265	LIGV(i,j)=cteurlig
266	POSLIGV(cteurlig,1)=i
267	POSLIGV(cteurlig,2)=j
268	ENDIF
269	END DO
270	END DO
271
272	! do i=1,cteurlig
273	! print *,'u',POSLIGU(i,1),POSLIGU(i,2),'v',POSLIGV(i,1),POSLIGV(i,2)
274	! enddo
275	! nblig est le nombre d'equation lineaire a resoudre dans mmat
276	if (NBLIG.ne.CTEURLIG) then
277	print *,"Erreur lors du second comptage"
278	stop
279	endif
280
281	942 continue
282	MMAT(:,:)=0
283	BDR(:,:)=0
284	BDRO(:)=0
285	!initialisation des tuij,suij... a 0.
286	TUMIMJ=0.; TUMIJ=0.; TUMIPJ=0.
287	TUIMJ=0. ; TUIJ=0. ; TUIPJ=0.
288	TUPIMJ=0.; TUPIJ=0.; TUPIPJ=0.
289	TVMIJ=0.;TVMIPJ=0.;TVIJ=0.;TVIPJ=0.
290
291	SVMIMJ=0.; SVMIJ=0.; SVMIPJ=0.
292	SVIMJ=0. ; SVIJ=0. ; SVIPJ=0.
293	SVPIMJ=0.; SVPIJ=0.; SVPIPJ=0.
294	SUIMJ=0.;SUIJ=0.;SUPIMJ=0.;SUPIJ=0.
295
296	!print*,'NBLIG =',NBLIG
297
298	! Initialisations et reperages
299	! ----------------------------
300	! IREORD=.FALSE.
301	! T pour faire le rearrangement de la matrice
302	! 0 pour aller le lire dans un fichier
303	! F pour traiter sans rearrangement
304
305	! On defini une valeur par defaut de kl et ku,
306	! 3_Boucle de determination des coefficients de la matrice
307	! ------------------------------------------------------
308	!!!!!!!!!! ATTENTION, le sens de l'enroulement pourrait etre modifié
309	DO I=1,NX
310	DO J=1,NY
311	!!!!!!!!!!!!!!!!!!!!!
312	!!!!!TEST POUR Uij!!!
313	!!!!!!!!!!!!!!!!!!!!!
314	IF ((I.eq.1.or.I.eq.NX).or.(J.eq.1.or.J.eq.NY)) THEN
315	TUIJ(I,J)=1
316	SVIJ(I,J)=1
317	METHOD(i,j,1)='bordx'
318	METHOD2(i,j,1)=-1
319	METHOD(i,j,2)='bordy'
320	METHOD2(i,j,2)=-1
321	ELSE
322	IF (.NOT.FLGZMX(I,J)) THEN !!Test sur FLGZMX
323	! points poses en x et n'etant pas sur le bord du domaine
324	TUIJ(I,J)=1
325	OPPOSX(i,j)=UXBAR(I,J)
326	if (okumat(i,j)) then
327	METHOD(i,j,1)='poseX'
328	METHOD2(i,j,1)=0
329	else
330	METHOD(i,j,1)='rienX'
331	METHOD2(i,j,1)=-1
332	endif
333	ELSE !!!Test sur FLGZMX, ici les points sont stream/shelf
334	!!on resout l'equation elliptique
335	!---------------------------------------------------
336
337	IF (ICE(i,j)==1) THEN!!!!!TEST SUR ICE, ici en glace
338
339	SELECT CASE (FRONT(i,j))
340	CASE (4) !le point a tous ses voisins => Calcul general
341	!
342	CALL U_GENERAL
343	METHOD(i,j,1)='4U_GEN'
344	METHOD2(i,j,1)=1
345	!
346	CASE (3) !le point a 3 voisins
347	IF (frontfacex(i,j)==0) THEN !!!!!!!!!!!!!!!!!!!!!!!front //Y
348	if (ICE(i-1,j+frontfacey(i,j))==1) then
349	CALL U_GENERAL
350	METHOD(i,j,1)='3U_GEN'
351	METHOD2(i,j,1)=1
352	else
353	CALL U_CISAILL(frontfacey(i,j))
354	METHOD(i,j,1)='3U_CIS'
355	METHOD2(i,j,1)=3
356	endif
357	ELSE !!!!!!!!!!!!!!!!!!!!!!!front //X
358	if (frontfacex(i,j)==-1) then !!front i-1 \| i (a droite)
359	CALL U_PRESSIO(-1)
360	METHOD(i,j,1)='3U_P-1'
361	METHOD2(i,j,1)=2
362	else
363	! CALL U_PRESSIO(+1)
364	CALL U_GENERAL
365	METHOD(i,j,1)='3U_GEN'
366	METHOD2(i,j,1)=1
367	endif
368	ENDIF
369	!
370	!
371	CASE (2) !le point a 2 voisins
372	IF (ISOLX(i,j)) THEN
373	!langue de glace non confinée avec bords //x \|\|
374	CALL U_TONGUE_LE(frontfacey(i,j))
375	METHOD(i,j,1)='2U_T_LE'
376	METHOD2(i,j,1)=6
377	ELSEIF (ISOLY(i,j)) THEN
378	!langue de glace non confinée avec bords //y =
379	CALL U_TONGUE
380	METHOD(i,j,1)='2U_T'
381	METHOD2(i,j,1)=5
382	ELSE ! ON A 2 FRONTS NON-OPPOSES
383	if (frontfacex(i,j)==-1) then
384	CALL U_COIN(-1)
385	METHOD(i,j,1)='2U_C-1'
386	METHOD2(i,j,1)=4
387	else
388	!!! CALL U_PRESSIO(+1)!pour Ui+1j!!!!!CALL U_COIN
389
390	IF((FRONT(i-1,j)==4).or.((FRONT(i-1,j)==3).and. &
391	(frontfacex(i-1,j)==-1)))then!!test pour Uij
392	CALL U_GENERAL
393	METHOD(i,j,1)='2U_GEN'
394	METHOD2(i,j,1)=1
395	ELSE
396	CALL U_CISAILL(frontfacey(i,j))
397	METHOD(i,j,1)='2U_CIS'
398	METHOD2(i,j,1)=3
399	ENDIF
400	endif
401	ENDIF
402	!
403	!
404	CASE (1)
405	IF (ISOLX(i,j)) THEN
406	!langue de glace non confinée avec bords //x \|\|
407	CALL U_TONGUE_LE(frontfacey(i,j))
408	METHOD(i,j,1)='1U_T_LE'
409	METHOD2(i,j,1)=6
410	ELSE
411	!langue de glace non confinée avec bords //y =
412	if (.not.ISOLY(i,j)) then!test
413	write(num_err,*) 'test sur u case=1 pb'!
414	STOP !test
415	endif !test
416	if (frontfacex(i,j)==1) then
417	CALL U_TONGUE
418	METHOD(i,j,1)='1U_T'
419	METHOD2(i,j,1)=5
420	else
421	CALL U_TONGUEFRONT(-1)
422	METHOD(i,j,1)='1U_TF-1'
423	METHOD2(i,j,1)=5
424	endif
425	ENDIF
426	CASE (0)
427	METHOD(i,j,1)='0U_'
428	METHOD2(i,j,1)=9
429	TUIJ(i,j)=1
430	OPPOSX(i,j)=5!
431	! PRINT*,'front =0','i=',i,'j=',j
432	CASE DEFAULT
433	write(num_err,*) 'ERREUR front est pas 0,1, 2, 3, 4'
434	STOP
435	END SELECT
436	ELSE!!!!!TEST SUR ICE, ij pas glace
437	IF (ISOLX(i-1,j)) THEN
438	CALL U_TONGUE_RI(frontfacey(i-1,j))
439	METHOD(i,j,1)='U_T_RI'
440	METHOD2(i,j,1)=6
441	ELSEIF (ISOLY(i-1,j)) THEN
442	CALL U_TONGUEFRONT(+1)
443	METHOD(i,j,1)='U_TF+1'
444	METHOD2(i,j,1)=5
445	ELSEIF (FRONT(i-1,j)==3) THEN
446	CALL U_PRESSIO(+1)
447	METHOD(i,j,1)='U_P+1'
448	METHOD2(i,j,1)=2
449	! TUIJ(i,j)=1
450	! opposx(i,j)=2
451	! print*, 'U_PRESSIO(+1) evite',i,j
452	ELSEIF (FRONT(i-1,j)==2) THEN
453	! if (ifail.gt.0.and.i==iii.and.j==jjj) then
454	! CALL U_PRESSIO(+1)
455	! METHOD(i,j,1)='U_P+1-'
456	! print*,'modif',i,j,METHOD(i,j,1)
457	if (ifail.gt.0.and.METHOD(i,j,1).eq.'U-chang') then
458	CALL U_PRESSIO(+1)
459	METHOD2(i,j,1)=10
460	! METHOD(i,j,1)='U_P+1-'
461	write(num_err,*)'modif',i,j,METHOD(i,j,1)
462	else
463	CALL U_COIN(+1)
464	METHOD2(i,j,1)=4
465	METHOD(i,j,1)='U_C+1'
466	endif
467	ELSE
468	METHOD(i,j,1)='U_ORFAN'
469	METHOD2(i,j,1)=999
470	TUIJ(i,j)=1
471	OPPOSX(i,j)=0!777666
472	! print*, 'point flgzmX sans glace', i,j
473	ENDIF
474
475	ENDIF!!!!!TEST SUR ICE, fin du test
476	ENDIF !!Test sur flgzmx
477
478	!!!!!!!!!!!!!!!!!!!!!
479	!!!!!TEST POUR Vij!!!
480	!!!!!!!!!!!!!!!!!!!!!
481	IF (.NOT.FLGZMY(I,J)) THEN !!Test sur FLGZMX/Y
482	SVIJ(I,J)=1
483	OPPOSY(i,j)=UYBAR(I,J)
484	!testvincent
485	if (okvmat(i,j)) then
486	METHOD(i,j,2)='poseY'
487	METHOD2(i,j,2)=0
488	else
489	METHOD(i,j,2)='rienY'
490	METHOD2(i,j,2)=-1
491	endif
492	ELSE !!!Test sur FLGZMY, ici les points sont stream/shelf
493	!! on resout l'equation elliptique
494	! ----------------------------------------------------
495	IF (ICE(i,j)==1) THEN !test sur ICE, ici point en glace
496	SELECT CASE (FRONT(i,j))
497	CASE (4) !le point a tous ses voisins => calcul general
498	!
499	CALL V_GENERAL
500	METHOD(i,j,2)='4V_GE'
501	METHOD2(i,j,2)=1
502	!
503	CASE (3) !le point a 3 voisins
504	!
505	IF (frontfacex(i,j)==0) THEN !!!!!!!!!!!!!!!!!!!!!!!front //Y
506
507	if (frontfacey(i,j)==-1) then !!front j-1 \| j en aval
508	CALL V_PRESSIO(-1)
509	METHOD(i,j,2)='3V_P-1'
510	METHOD2(i,j,2)=2
511	else !!front j \| j+1 en amont
512	! CALL V_PRESSIO(+1)
513	CALL V_GENERAL
514	METHOD(i,j,2)='3V_GE'
515	METHOD2(i,j,2)=1
516	endif
517	ELSE !!!!!!!!!!!!!!!!!!!!!!!front //X
518	!test pour le calcul de Vij
519	if(ICE(i+frontfacex(i,j),j-1)==1) then
520	CALL V_GENERAL
521	METHOD(i,j,2)='3V_GEN'
522	METHOD2(i,j,2)=1
523	else
524	CALL V_CISAILL(frontfacex(i,j))
525	METHOD(i,j,2)='3V_CIS'
526	METHOD2(i,j,2)=3
527	endif
528	ENDIF
529	!
530	!
531	CASE (2)
532	!
533	IF (ISOLX(i,j)) THEN
534	!langue de glace non confinée avec bords //x \|\|
535	CALL V_TONGUE
536	METHOD(i,j,2)='2V_T'
537	METHOD2(i,j,2)=5
538	ELSEIF (ISOLY(i,j)) THEN
539	!langue de glace non confinée avec bords //y =
540	CALL V_TONGUE_DO(frontfacex(i,j))
541	METHOD(i,j,2)='2V_T_DO'
542	METHOD2(i,j,2)=6
543	ELSE ! ON A 2 FRONTS NON-OPPOSES
544
545	if (frontfacey(i,j)==-1) then
546	CALL V_COIN(-1)
547	METHOD(i,j,2)='2V_C-1'
548	METHOD2(i,j,2)=4
549	else
550	!!!!!!!!!!!!! CALL V_PRESSIO(+1)!!!!!!!!!!!!!!!!CALL V_COIN
551	!TESTS POUR VOIR DE QUEL TYPE EST VIJ
552
553	IF((FRONT(i,j-1)==4).or.((FRONT(i,j-1)==3).and. &
554	(frontfacey(i,j-1)==-1)))then
555	CALL V_GENERAL
556	METHOD(i,j,2)='2V_GE'
557	METHOD2(i,j,2)=1
558	ELSE
559	CALL V_CISAILL(frontfacex(i,j))
560	METHOD(i,j,2)='2V_CIS'
561	METHOD2(i,j,2)=3
562	ENDIF
563	endif
564	ENDIF
565	!
566	CASE (1)
567	IF (ISOLX(i,j)) THEN
568	!!!!!!!!!!!!!!!!!!! if (frontfacex(i,j).eq.1) then
569	if (frontfacey(i,j).eq.1) then
570	CALL V_TONGUE
571	METHOD(i,j,2)='1V_T'
572	METHOD2(i,j,2)=5
573	else
574	CALL V_TONGUEFRONT(-1)
575	METHOD(i,j,2)='1V_TF-1'
576	METHOD2(i,j,2)=5
577	endif
578	ELSEIF (ISOLY(i,j)) THEN
579	CALL V_TONGUE_DO(frontfacex(i,j))
580	METHOD(i,j,2)='1V_T_DO'
581	METHOD2(i,j,2)=6
582	ELSE
583	write(num_err,*) 'ERREUR test sur case 1 ds remplimat'
584	STOP
585	ENDIF
586	CASE (0)
587	METHOD(i,j,2)='0V_'
588	METHOD2(i,j,2)=9
589	SVIJ(i,j)=1
590	OPPOSY(i,j)=5!
591	CASE DEFAULT
592	write(num_err,*) 'ERREUR front est pas 0,1, 2, 3, 4'
593	STOP
594	END SELECT
595	ELSE !!!!test sur ice(i,j) point non en glace
596	if (isolx(i,j-1)) then
597	CALL V_TONGUEFRONT(+1)
598	METHOD(i,j,2)='V_TF+1'
599	METHOD2(i,j,2)=5
600	elseif (isoly(i,j-1)) then
601	CALL V_TONGUE_UP(frontfacex(i,j-1))
602	METHOD(i,j,2)='V_T_UP'
603	METHOD2(i,j,2)=6
604	elseif (front(i,j-1)==3) then
605	if (ifail.gt.0.and.METHOD(i,j,2).eq.'V-chang') then
606	CALL U_COIN(+1)
607	METHOD2(i,j,1)=10
608	write(num_err,*) 'modif',i,j,METHOD(i,j,1)
609	else
610	CALL V_PRESSIO(+1)
611	METHOD(i,j,2)='V_P+1'
612	METHOD2(i,j,2)=2
613	endif
614	elseif (front(i,j-1)==2) then
615	CALL V_COIN(+1)
616	METHOD(i,j,2)='V_C+1'
617	METHOD2(i,j,2)=4
618	else
619	METHOD(i,j,2)='V_ORFAN'
620	METHOD2(i,j,2)=995
621	SVIJ(i,j)=1
622	OPPOSY(i,j)=0!888666
623	! print*, 'point flgzmY sans ice',i,j
624	endif
625
626	ENDIF !!test sur ice(i,j), fin du test
627	ENDIF !!Test sur flgzmy
628	ENDIF
629	ENDDO
630	ENDDO
631
632	write(num_err,*) 'remplimat he'
633	! 4_Boucle de remplissage de la matrice.
634	! -------------------------------------
635	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
636	! on determine la position de chaque vitesse Uij et Vij
637	!!! dec 2003 pour l'instant on considere tous les points
638
639	! DO J=1,NY
640	! DO I=1,NX
641	! LIGU(i,j)=CTEURLIG
642	! LIGV(i,j)=CTEURLIG+1
643	! cteurlig=cteurlig+2
644	! ENDDO
645	! ENDDO
646	! nblig=cteurlig-1
647	!stop
648	iinf=51
649	isup=53
650	jinf=53
651	jsup=55
652	CTEURLIG = 0
653	aaa=0.
654	! print*,'remplissage de mmat'
655	! print*,('shape=',shape(mmat))
656	! print*,('size =',size(mmat))
657	! print*,'remplissage de bdr'
658	! print*,('shape=',shape(bdr))
659	! print*,('size =',size(bdr))
660
661	DO I=1,NX
662	DO J=1,NY
663
664	test_oku : IF (OKUMAT(I,J)) THEN
665	CTEURLIG = CTEURLIG + 1
666	! print*,i,j,'okumat=',OKUMAT(I,J),'CTEUR=',CTEURLIG,'LIGU',LIGU(I,J)
667
668	test_shelfx : IF ((I.eq.1.or.I.eq.NX).or.(J.eq.1.or.J.eq.NY) &
669	.or.(.not.FLGZMX(I,J))) THEN
670	mmat(kdc2,ligu(i,j)) = 1.0
671	! BDR(ligu(i,j),1) =UXBAR(I,J)
672	BDR(ligu(i,j),1) =OPPOSX(I,J)
673
674
675	if (OPPOSX(I,J).ne.UXBAR(I,J)) then
676	write(num_err,*) 'pb OPPOSX pose',i,j
677	write(num_err,*) METHOD(i,j,1),OPPOSX(I,J),UXBAR(I,J)
678
679
680
681
682
683
684	endif
685	ELSE
686
687	! Classement de la matrice pour Uij
688	! ---------------------------------
689
690	SCAL = TUIJ(i,j)
691
692	mmat(kdc2,ligu(i,j)) = 1.0
693	iimat=kdc2+ligu(i,j)
694	!print*,'uij mmat(kdc2,ligu)=mmat(',kdc2,ligu(i,j),')'
695	BDR(ligu(i,j),1) = OPPOSX(i,j)/SCAL
696
697	opposx(i,j) = OPPOSX(i,j)/SCAL
698
699	if (LIGU(I-1,J-1).gt.0) then
700	SNR = LIGU(I-1,J-1) ! les Ui-1
701	!print*, 'mmat(',iimat-snr,snr,') TUMIMJ'
702	!print*,'iimat=',iimat,'snr=',snr
703	mmat(iimat-snr,snr) = TUMIMJ(i,j)/SCAL
704	endif
705	if (LIGU(I-1,J).gt.0) then
706	SNR = LIGU(I-1,J)
707	mmat(iimat-snr,snr) = TUMIJ(i,j)/SCAL
708	endif
709	if (LIGU(I-1,J+1).gt.0) then
710	SNR = LIGU(I-1,J+1)
711
712	mmat(iimat-snr,snr) = TUMIPJ(i,j)/SCAL
713	endif
714	if (LIGU(I,J-1).gt.0) then
715	SNR = LIGU(I,J-1) ! les Ui
716	mmat(iimat-snr,snr) = TUIMJ(i,j)/SCAL
717	endif
718	if (LIGU(I,J+1).gt.0) then
719	SNR = LIGU(I,J+1)
720	mmat(iimat-snr,snr) = TUIPJ(i,j)/SCAL
721	endif
722	if (LIGU(I+1,J-1).gt.0) then
723	SNR = LIGU(I+1,J-1) ! les Ui+1
724	mmat(iimat-snr,snr) = TUPIMJ(i,j)/SCAL
725	endif
726	if (LIGU(I+1,J).gt.0) then
727	SNR = LIGU(I+1,J)
728
729	mmat(iimat-snr,snr) = TUPIJ(i,j)/SCAL
730	endif
731	if (LIGU(I+1,J+1).gt.0) then
732	SNR = LIGU(I+1,J+1)
733	mmat(iimat-snr,snr) = TUPIPJ(i,j)/SCAL
734
735	endif
736	if (LIGV(I-1,J).gt.0) then
737	SNR = LIGV(I-1,J) ! les Vi-1
738	mmat(iimat-snr,snr) = TVMIJ(i,j)/SCAL
739
740	endif
741	if (LIGV(I-1,J+1).gt.0) then
742	SNR = LIGV(I-1,J+1)
743	mmat(iimat-snr,snr) = TVMIPJ(i,j)/SCAL
744	endif
745	if (LIGV(I,J).gt.0) then
746	SNR = LIGV(I,J) ! les Vi
747	mmat(iimat-snr,snr) = TVIJ(i,j)/SCAL
748	endif
749	if (LIGV(I,J+1).gt.0) then
750	SNR = LIGV(I,J+1)
751	!print*, 'mmat(',iimat-snr,snr,') TVIPJ'
752	!print*,'iimat=',iimat,'snr=',snr
753	mmat(iimat-snr,snr) = TVIPJ(i,j)/SCAL
754	endif
755	ENDIF test_shelfx
756	ENDIF test_oku
757
758
759	test_okv : IF (OKVMAT(I,J)) THEN
760	CTEURLIG = CTEURLIG + 1
761	! print*,i,j,'okvmat=',OKUMAT(I,J),'CTEUR=',CTEURLIG,'LIGV',LIGV(I,J)
762
763	test_shelfy : IF ((I.eq.1.or.I.eq.NX).or.(J.eq.1.or.J.eq.NY) &
764	.or.(.not.FLGZMY(I,J))) THEN
765	mmat(kdc2,ligv(i,j)) = 1.0
766	! BDR(ligv(i,j),1) =UYBAR(I,J)
767	BDR(ligv(i,j),1) =OPPOSY(I,J)
768	ELSE !test_shelfy
769	! Classement de la matrice pour Vij
770	! ---------------------------------
771
772	SCAL = SVIJ(i,j)
773
774	iimat=kdc2+ligv(i,j)
775	!print*,'vij mmat(kdc2,ligv)=mmat(',kdc2,ligv(i,j),')'
776	mmat(kdc2,ligv(i,j)) = 1.0
777
778	BDR(ligv(i,j),1) = OPPOSY(i,j)/SCAL
779	if ((i.ge.67).and.(i.le.71)) then
780	write(6,*)'i,j,opposy,scal',i,j,OPPOSY(i,j),scal,OPPOSY(i,j)/scal,method(i,j,2)
781	end if
782
783	opposy(i,j)=opposy(i,j)/scal
784
785	if (LIGV(I-1,J-1).gt.0) then
786	SNR = LIGV(I-1,J-1) ! les Vi-1
787	mmat(iimat-snr,snr) = SVMIMJ(i,j)/SCAL
788	endif
789
790	if (LIGV(I-1,J).gt.0) then
791	SNR = LIGV(I-1,J)
792	mmat(iimat-snr,snr) = SVMIJ(i,j)/SCAL
793	endif
794
795	if (LIGV(I-1,J+1).gt.0) then
796	SNR = LIGV(I-1,J+1)
797	!print*, 'mmat(',iimat-snr,snr,') SVMIPJ'
798	!print*,'iimat=',iimat,'snr=',snr
799	mmat(iimat-snr,snr) = SVMIPJ(i,j)/SCAL
800	endif
801
802	if (LIGV(I,J-1).gt.0) then
803	SNR = LIGV(I,J-1) ! les Vi
804	mmat(iimat-snr,snr) = SVIMJ(i,j)/SCAL
805	endif
806
807	if (LIGV(I,J+1).gt.0) then
808	SNR = LIGV(I,J+1)
809	mmat(iimat-snr,snr) = SVIPJ(i,j)/SCAL
810	endif
811
812	if (LIGV(I+1,J-1).gt.0) then
813	SNR = LIGV(I+1,J-1) ! les Vi+1
814	mmat(iimat-snr,snr) = SVPIMJ(i,j)/SCAL
815	endif
816
817	if (LIGV(I+1,J).gt.0) then
818	SNR = LIGV(I+1,J)
819	mmat(iimat-snr,snr) = SVPIJ(i,j)/SCAL
820	endif
821
822	if (LIGV(I+1,J+1).gt.0) then
823	SNR = LIGV(I+1,J+1)
824
825	mmat(iimat-snr,snr) = SVPIPJ(i,j)/SCAL
826	endif
827
828	if (LIGU(I,J-1).gt.0) then
829	SNR = LIGU(I,J-1) ! les Ui
830	mmat(iimat-snr,snr) = SUIMJ(i,j)/SCAL
831	endif
832	if (LIGU(I,J).gt.0) then
833	SNR = LIGU(I,J)
834	mmat(iimat-snr,snr) = SUIJ(i,j)/SCAL
835	endif
836
837	if (LIGU(I+1,J-1).gt.0) then
838	SNR = LIGU(I+1,J-1) !les Ui+1
839	mmat(iimat-snr,snr) = SUPIMJ(i,j)/SCAL
840	endif
841
842	if (LIGU(I+1,J).gt.0) then
843	SNR = LIGU(I+1,J)
844	mmat(iimat-snr,snr) = SUPIJ(i,j)/SCAL
845	endif
846	ENDIF test_shelfy
847	ENDIF test_okv
848
849	!print*,('classement des V fini')
850
851	ENDDO
852	ENDDO
853	!print*,'CTEURLIG',CTEURLIG
854	!print*,'NBLIG =',NBLIG
855	NBLIG = CTEURLIG - 1
856
857	!print*,'fin de boucle de remplissage'
858	! stop
859	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
860	! 5_Resolution du systeme lineaire
861	! ------------------------------
862	ifail=0
863	conjgr=.false.
864
865	if (CONJGR) then !!!!!!!!!!!!debut test sur CONJGR
866	DO K=1,NBLIG
867	BDRO(K) = BDR(K,1)
868	END DO
869	! nouveau conjgrad (utilisant blas)
870	! -------------------------
871	! ----------------a enlever apres test-----------------------
872	CONJGR = .true.
873	DO I=1,NX
874	DO J=1,NY
875	! DO I=ID1,ID2
876	! DO J=JD1,JD2
877	IF (OKUMAT(I,J)) THEN
878	SAV(LIGU(I,J))=UXPREC(I,J)
879	END IF
880	IF (OKVMAT(I,J)) THEN
881	SAV(LIGV(I,J))=UYPREC(I,J)
882	END IF
883	END DO
884	END DO
885	! -----------------fin de la partie a enlever---------------
886	! write(6,*)'avant conjgrad time=',time,' ldab=',ldab
887	! call conjgrad(nblig,kl,ku,mmat,ldab,BDR,SAV,rsq,iter,ifail)
888	! write(6,*)'time',time,' numdom',numdom,' iter',iter
889
890	! CALL CONJGRAD(BDRO,NBLIG,SAV,RSQ,ifail,ITER)
891	if (iter.gt.1) write(77,*) 'time=',time,' iter=',iter
892	IF (ifail.EQ.1) THEN
893	write(,) 'gradients conjugues insuffisants'
894	CONJGR=.FALSE.
895	TIMECG=TIME
896	goto 911
897
898	ENDIF
899
900	ELSE !!!!!!!!!!!!continue test sur CONJGR
901
902	! write(6,*)
903	! write(6,*) 'avant entree sgbsv, domaine',numdom
904	! write(6,*) 'nblig',nblig,' kl',kl,' ku',ku,' nrhs',nrhs
905	! write(6,*) 'ldab',ldab,' ldb',ldb
906
907	knonul=0
908	emmax=0.
909
910	!----pour mettre mmat=identite----!
911	! do i=1,ldbmax
912	! do j=1,nptmax
913	! if (i.eq.ku+kl+1) then
914	! mmat(i,j)=1.
915	! else
916	! mmat(i,j)=0.
917	! endif
918	! end do
919	! end do
920	!---------------------------------!
921
922	do i=1,nblig
923	bdro(i)=bdr(i,1)
924	end do
925
926	! pour verifier qu'on a pas oublie de point
927	do i=1,nblig
928	if (mmat(kl+ku+1,i).eq.0.) then
929	print*,(mmat(kl+ku+1,i))
930	print*,'kl+ku+1',kl+ku+1
931	write(6,) '****** ATTENTION point oublie numero' &
932	,i,'****************************'
933	endif
934	end do
935	!----chercher le maximum de la matrice-----------------!
936	! emmax=0.
937	! do i=1,nblig
938	! do j=1,ldab
939	! emmax=max(abs(mmat(j,i)),emmax)!-!
940	! end do
941	! end do
942	! write(6,*) emmax
943	!------------------------------------------------------!
944	!open(UNIT=994,file='indexes')
945	! do i=1,nx
946	! do j=1,ny
947	! write(994,*) 'i ,j ,LIGU,LIGV'
948	! write(994,*) i,j,LIGU(i,j),LIGV(i,j)
949	! enddo
950	! enddo
951	! close (994)
952	!open(UNIT=995,file='bdr')
953	! do cteurlig=1,nlig
954	! write(995,*) bdr(cteurlig),cteurlig
955	! enddo
956	! close (995)
957	!open(UNIT=996,file='mmat_columns_before')
958	!write(996,*)'bdr(u(10,10))=',bdr(LIGU(10,10))
959	!write(996,*)'bdr(v(10,10))=',bdr(LIGV(10,10))
960	!do i=1,2*ku+ku+1
961	! write(996,*) mmat(i,LIGU(10,10)),mmat(i,LIGV(10,10))
962	!enddo
963	!close (996)
964	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
965	! write(6,*)'avant sgbsv time=',time,' ldab=',ldab
966	! write(6,*)'nblig',nblig,'ku=',ku,'kl=',kl,'nrhs',nrhs
967	! write(6,*)'ipiv ?','ldb ',ldb,'ifail',ifail
968	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
969	! print*, 'min value, location',minval(mmat),minloc(mmat)
970	! print*, 'max value, location',maxval(mmat),maxloc(mmat)
971	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
972
973
974	! sorties des Tuij
975	debug_3D(:,:,26)=TUIJ(:,:)
976
977
978
979	! write(6,*)'avant sgbsv', nblig,kl,ku,nrhs,ldbmax
980
981	! quelques sorties netcdf
982
983	where (okumat(:,:))
984	debug_3D(:,:,39)=1
985	elsewhere
986	debug_3D(:,:,39)=0
987	end where
988
989	where (okvmat(:,:))
990	debug_3D(:,:,40)=1
991	elsewhere
992	debug_3D(:,:,40)=0
993	end where
994
995
996	debug_3d(:,:,35) = opposx(:,:)
997	debug_3d(:,:,36) = opposy(:,:)
998
999
1000
1001
1002	CALL SGBSV(nblig,kl,ku,nrhs,mmat,ldbmax,ipiv,bdr,ldb,ifail)
1003	! write(6,*)'apres sgbsv'
1004	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1005	! open(UNIT=998,file='sgbsv_vit_u.dat')
1006	! open(UNIT=999,file='sgbsv_vit_v.dat')
1007	! do i=1,nblig
1008	! if (bdro(i).ne.bdr(i,1)) then
1009	!write(6,*) i, bdro(i)-bdr(i,1)
1010	! endif
1011	! end do
1012
1013	! call DEBUGAGE
1014
1015	! write(num_ic_by,*) mmat
1016	! write(6,*) 'fin de l ecriture de mmat'
1017
1018
1019	IF (ifail.NE.0) THEN
1020	nb_err=nb_err+1
1021	WRITE(,) ' ======================================'
1022	WRITE(,) ' ERREUR DANS L''ELIMINATION GAUSSIENNE '
1023	WRITE(,) ' ERREUR N ',ifail
1024	WRITE(,) ' ======================================'
1025
1026	if (posligu(ifail,1).gt.0) then
1027	iii=posligu(ifail,1)
1028	jjj=posligu(ifail,2)
1029	print*,'Ux',iii,jjj,method(iii,jjj,1)
1030	METHOD(iii,jjj,1)='U-chang'
1031	elseif (posligv(ifail,1).gt.0) then
1032	iii=posligv(ifail,1)
1033	jjj=posligv(ifail,2)
1034	print*,'Vy',iii,jjj,method(iii,jjj,2)
1035	METHOD(iii,jjj,2)='V-chang'
1036	endif
1037	if (geoplace.eq.'anteis1'.and.H(71,71).lt.10.) then
1038	print*,time,ifail
1039	print*,'stop apres glace disparut a Pole sud h=',h(71,71)
1040	stop
1041	endif
1042	call DEBUGAGE
1043	!!! if (nb_err.lt.2) then
1044	!!! goto 942
1045	!!! else
1046	print*,'=apres',nb_err,'passage infructueux : on stop'
1047	! pause
1048	return
1049	! STOP
1050	!!! endif
1051	END IF
1052
1053	! call DEBUGAGE
1054	DO I=1,NX
1055	DO J=1,NY
1056
1057	! if ((I.ge.ID1).and.(I.le.ID2).and. &
1058	! (J.ge.JD1).and.(J.le.JD2).and.(OKUMAT(I,J))) then
1059	if (OKUMAT(I,J)) then
1060	UXNEW(I,J)=BDR(LIGU(I,J),1)
1061	SAV(LIGU(I,J))=UXNEW(I,J)
1062	else
1063	UXNEW(I,J)=UXPREC(I,J)
1064	endif
1065
1066	! if ((I.ge.ID1).and.(I.le.ID2).and. &
1067	! (J.ge.JD1).and.(J.le.JD2).and.(OKVMAT(I,J))) then
1068	if (OKVMAT(I,J)) then
1069	UYNEW(I,J)=BDR(LIGV(I,J),1)
1070	SAV(LIGV(I,J))=UYNEW(I,J)
1071	else
1072	UYNEW(I,J)=UYPREC(I,J)
1073	endif
1074	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1075	!open(UNIT=998,file='sgbsv_vit_u.dat')
1076	!open(UNIT=999,file='sgbsv_vit_v.dat')
1077	! write (998,*) i,j,UXNEW(I,J)
1078	! write (999,*) i,j,UYNEW(I,J)
1079	if (ICE(i,j)==0) then
1080	if (.not.flgzmx(i,j)) then
1081	UXNEW(I,J)=3.3333333
1082	endif
1083	if (.not.flgzmy(i,j)) then
1084	UYNEW(I,J)=3.3333333
1085	endif
1086	endif
1087	!close (998)
1088	!close (999)
1089	!!!!!!!!!!!!!test vincent : limiter vitesses apres sgbsv
1090	UXNEW(I,J)=max(-3900.,UXNEW(I,J))
1091	UXNEW(I,J)=min( 3900.,UXNEW(I,J))
1092	UYNEW(I,J)=max(-3900.,UYNEW(I,J))
1093	UYNEW(I,J)=min( 3900.,UYNEW(I,J))
1094	END DO
1095	END DO
1096
1097	where ((Hmx(:,:).le.1.).and.(flgzmx(:,:)))
1098	uxnew(:,:)=0.
1099	end where
1100
1101	where ((Hmy(:,:).le.1.).and.(flgzmy(:,:)))
1102	uynew(:,:)=0.
1103	end where
1104
1105
1106	debug_3D(:,:,27)=uxnew(:,:)
1107
1108
1109
1110
1111	! print*,'???',UXNEW(57,40),UXNEW(58,40),UXNEW(59,40)
1112	! pause
1113	! close (998)
1114	!!!!!!! close (999)
1115	ENDIF !!!!!!!!!!!!fin test sur CONJGR
1116
1117	! impression des vitesses
1118
1119	! open (57,file='Dom'//domname//domtot2(1:jn))
1120	! write(57,*) 'time=',time
1121	! write(57,*)'i j vitesse old-vitesse delta-vitesse '
1122	! do i=2,nx-1
1123	! do j=2,ny-1
1124	! UNW = SQRT(((UXnew(I,J)+UXnew(I+1,J))/2)**2
1125	! & + ((UYnew(I,J)+UYnew(I,J+1))/2.)**2)
1126	! UPR = SQRT(((UXPREC(I,J)+UXPREC(I+1,J))/2)**2
1127	! & + ((UYPREC(I,J)+UYPREC(I,J+1))/2.)**2)
1128
1129	! write(57,958) i,j,unw, upr,unw-upr
1130	! end do
1131	! end do
1132	! close(57)
1133	!958 format(i3,1x,i3,2x,3(e10.4,2x))
1134	! print*,'fin de remplimat'
1135
1136	911 continue
1137
1138	close(num_err)
1139
1140
1141	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1142
1143	! 6_routines pour referencer les coefficients utiles pour remplir la matrice
1144	! ------------------------------------------------------------------------
1145	CONTAINS ! Les 10 procedures suivantes sont internes a remplidom
1146	! U_* V_* determinent les procedures de calcul du shelf
1147	! Cas general U_GENERAL V_GENERAL
1148	! face a un front action (normale au front) de la pression U_PERSSION V_PERSSION
1149	! en bordure du font condition de non cisaillzment U_CISAILL V_CISAILL
1150	! langue non confinee ISOLXBACKW ISOLYBACKW calcul des vitesses en backward
1151	! languee non confinee ISOLXFOREW ISOLYFOREW calcul des
1152	! vitesses en foreward
1153	! LARGEUR_DE_BANDE determine les ku et kl minimums
1154	! DEBUG ecrit une serie de fichiers TABLE_qch qui donne les coeff de
1155	! mmat, la methode utilisees
1156	SUBROUTINE U_GENERAL
1157	!Rempli la matice pour le calcul de Uij
1158	IMPLICIT NONE
1159
1160
1161	! Calcul des coefficients de l'equation en U:
1162	! -------------------------------------------
1163	beta=FROTMX(I,J)dxdx
1164	! Terme en u(i,j):
1165	! _______________
1166	TUIJ(I,J) = -4.PVI(I,J) - 4.PVI(I-1,J) &
1167	- PVM(i,j+1)-PVM(I,J)-beta
1168	! Terme en u(i-1,j):
1169	! _________________
1170	TUMIJ(I,J) = 4.*PVI(I-1,J)
1171
1172	! Terme en u(i+1,j):
1173	! _________________
1174	TUPIJ(I,J) = 4.*PVI(I,J)
1175
1176	! Terme en u(i,j-1):
1177	! _________________
1178	TUIMJ(I,J) = PVM(I,J)
1179
1180	! Terme en u(i,j+1):
1181	! _________________
1182	TUIPJ(I,J) = PVM(I,J+1)
1183
1184	! Terme en v(i,j):
1185	! _______________
1186	TVIJ(I,J) = -2.*PVI(I,J)-PVM(I,J)
1187
1188	! Terme en v(i-1,j):
1189	! _________________
1190	TVMIJ(I,J) = 2.*PVI(I-1,J)+PVM(I,J)
1191
1192	! Terme en v(i-1,j+1):
1193	! ___________________
1194	TVMIPJ(I,J) = -2.*PVI(I-1,J)-PVM(I,J+1)
1195
1196	! Terme en v(i,j+1):
1197	! _________________
1198	TVIPJ(I,J) = 2.*PVI(I,J)+PVM(I,J+1)
1199
1200	moteur= ROGHMX(I,J)* &
1201	(S(I,J)-S(I-1,J))/dx
1202	moteur=min(moteur,moteurmax)
1203	moteur=max(moteur,-moteurmax)
1204
1205	OPPOSX(I,J) = moteurdxdx
1206	!print*,'u_gen OPPOS',OPPOSX(I,J),i,j
1207	!print,'motmax=',moteurmax,'motmaxdx2=',moteurmaxdxdx
1208
1209	END SUBROUTINE U_GENERAL
1210	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1211	SUBROUTINE V_GENERAL
1212	!Rempli la matice pour le calcul de Vij
1213	IMPLICIT NONE
1214
1215	! print*, 'V_GENERAL',i,j
1216
1217	! Calcul des coefficients de l'equation en V:
1218	! -------------------------------------------
1219	beta=FROTMY(I,J)dxdx
1220
1221	! Terme en v(i,j):
1222	! ________________
1223	SVIJ(I,J) = -4.PVI(I,J)-4.PVI(I,J-1) &
1224	- PVM(I+1,J)-PVM(I,J)-beta
1225
1226	! Terme en v(i,j-1):
1227	! __________________
1228	SVIMJ(I,J) = 4.*PVI(I,J-1)
1229
1230	! Terme en v(i,j+1):
1231	! __________________
1232	SVIPJ(I,J) = 4.*PVI(I,J)
1233
1234	! Terme en v(i-1,j):
1235	! __________________
1236	SVMIJ(I,J) = PVM(I,J)
1237
1238	! Terme en v(i+1,j):
1239	! __________________
1240	SVPIJ(I,J) = PVM(I+1,J)
1241
1242	! Terme en u(i,j):
1243	! ________________
1244	SUIJ(I,J) = -2.*PVI(I,J)-PVM(I,J)
1245
1246	! Terme en u(i,j-1):
1247	! __________________
1248	SUIMJ(I,J) = 2.*PVI(I,J-1)+PVM(I,J)
1249
1250
1251	! Terme en u(i+1,j-1):
1252	! ___________________
1253	SUPIMJ(I,J) = -2.*PVI(I,J-1)-PVM(I+1,J)
1254
1255
1256	! Terme en u(i+1,j):
1257	! ___________________
1258	SUPIJ(I,J) = 2.*PVI(I,J)+PVM(I+1,J)
1259
1260
1261	moteur= ROGHMY(I,J)* &
1262	(S(I,J)-S(I,J-1))/dx
1263	moteur=min(moteur,moteurmax)
1264	moteur=max(moteur,-moteurmax)
1265
1266	OPPOSY(i,j) = moteurdxdx
1267
1268
1269	END SUBROUTINE V_GENERAL
1270	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1271
1272
1273	SUBROUTINE U_PRESSIO(DIR)
1274	! sert quand front=3 et frontfacex=+-1
1275	!Rempli la matice pour le calcul de Uij
1276	IMPLICIT NONE
1277
1278	INTEGER DIR!Direction du front vers l'aval+1,vers l'amont-1
1279	! Equivaut frontfacey
1280
1281	if (DIR.EQ.1) then !front a droite (case vide a droite de Uij)
1282	beta=FROTMX(I,J)*dx
1283	beta=min(beta,0.34pvi(I-1,J))
1284	TUIJ(I,J) = 4*PVI(I-1,J) - beta
1285	TUMIJ(I,J) =-4*PVI(I-1,J)
1286	TVMIPJ(I,J) = 2*PVI(I-1,J)
1287	TVMIJ(I,J) =-2*PVI(I-1,J)
1288	!!!! OPPOSX(I,J) = ROG(1.-RO/ROW)(H(I-1,J)*2)dx/2.
1289	OPPOSX(I,J) = ROG(H(I-1,J)2)dx/2. - ROWG((min(slv(i-1,j)-B(i-1,j),0.))2)dx/2.
1290	else !front a gauche (case vide a gauche de Uij)
1291	beta=FROTMX(I,J)*dx
1292	beta=min(beta,0.34pvi(I,J))
1293	TUIJ(I,J) =-4*PVI(I,J) + beta
1294	TUPIJ(I,J) = 4*PVI(I,J)
1295	TVIPJ(I,J) = 2*PVI(I,J)
1296	TVIJ(I,J) =-2*PVI(I,J)
1297	!!!! OPPOSX(I,J) = ROG(1.-RO/ROW)(H(I,J)*2)dx/2.
1298	OPPOSX(I,J) = ROG(H(I,J)2)dx/2. - ROWG((min(slv(i,j)-B(i,j),0.))2)dx/2.
1299	endif
1300
1301	END SUBROUTINE U_PRESSIO
1302	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1303	SUBROUTINE V_PRESSIO(DIR)
1304	! sert quand front=3 et frontfacey=+-1
1305	!Rempli la matice pour le calcul de Vij
1306	IMPLICIT NONE
1307
1308	INTEGER DIR!Direction du front vers l'aval+1,vers l'amont-1
1309	! Equivaut frontfacey
1310
1311	if (DIR.EQ.1) then!front en haut (case vide EN HAUT de Vij)
1312
1313	beta=FROTMY(I,J)*dx
1314	beta=min(beta,0.34pvi(I,J-1))
1315	SVIJ(I,J) = 4*pvi(I,J-1) - beta
1316	SVIMJ(I,J) =-4*pvi(I,J-1)
1317	SUPIMJ(I,J) = 2*pvi(I,J-1)
1318	SUIMJ(I,J) =-2*pvi(I,J-1)
1319	OPPOSY(I,J) = ROG(H(I,J-1)2)dx/2. - ROWG((min(slv(i,j-1)-B(i,j-1),0.))2)dx/2.
1320	else!front en bas (case vide en bas de Vij)
1321	! print*, 'V_PRESSIO(-1)',i,j
1322	beta=FROTMY(I,J)*dx
1323	beta=min(beta,0.34pvi(I,J))
1324	SVIJ(I,J) =-4*pvi(I,J) + beta
1325	SVIPJ(I,J) = 4*pvi(I,J)
1326	SUPIJ(I,J) = 2*pvi(I,J)
1327	SUIJ(I,J) =-2*pvi(I,J)
1328	!!!! OPPOSY(I,J) = ROG(1.-RO/ROW)(H(I,J)*2)dx/2.
1329	OPPOSY(I,J) = ROG(H(I,J)2)dx/2. - ROWG((min(slv(i,j)-B(i,j),0.))2)dx/2.
1330
1331	end if
1332
1333	END SUBROUTINE V_PRESSIO
1334	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1335	SUBROUTINE U_COIN(CAS)
1336	! sert quand front=2 et frontfacex=+-1,frontfacey=+-1
1337	!Rempli la matice pour le calcul de Uij
1338	IMPLICIT NONE
1339
1340	integer CAS !valeur de frontfacex :si le front est a-droite ou
1341	!a-gauche du point en glace
1342	if (cas==1) then
1343	beta=FROTMX(i,j)*dx
1344	beta=min(beta,0.3*PVI(i-1,j))
1345	! TUIJ(I,J) = 6*PVI(i-1,j) - beta
1346	! TUMIJ(I,J) =-6*PVI(i-1,j)
1347	TUIJ(I,J) = PVI(i-1,j) - beta
1348	TUMIJ(I,J) =-PVI(i-1,j)
1349	!!!! OPPOSX(I,J) = ROG(1.-RO/ROW)(H(I-1,J)*2)dx/2.
1350	OPPOSX(i,j)= ROG(H(I-1,J)2)dx/2. - ROWG((min(slv(i-1,j)-B(i-1,j),0.))2)dx/2.
1351
1352	else
1353	beta=FROTMX(i,j)*dx
1354	beta=min(beta,0.3*PVI(i,j))
1355	! TUIJ(I,J) =-6*PVI(i,j) + beta
1356	! TUPIJ(I,J) = 6*PVI(i,j)
1357	TUIJ(I,J) =-PVI(i,j) + beta
1358	TUPIJ(I,J) = PVI(i,j)
1359	!!!! OPPOSX(I,J) = ROG(1.-RO/ROW)(H(I,J)*2)dx/2.
1360	OPPOSX(i,j)= ROG(H(I,J)2)dx/2.- ROWG((min(slv(i,j)-B(i,j),0.))2)dx/2.
1361
1362	endif
1363	END SUBROUTINE U_COIN
1364	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1365	SUBROUTINE V_COIN(CAS)
1366	! sert quand front=2 et frontfacex=+-1,frontfacey=+-1
1367	!Rempli la matice pour le calcul de Vij
1368	IMPLICIT NONE
1369	integer CAS !valeur de frontfacey : si le front est au-dessus ou
1370	!au-dessous du point en glace
1371	if (cas==1) then
1372	beta=FROTMY(I,J)*dx
1373	beta=min(beta,0.3*PVI(i,j-1))
1374	! SVIJ(I,J) = 6*PVI(i,j-1) - beta
1375	! SVIMJ(I,J) =-6*PVI(i,j-1)
1376	SVIJ(I,J) = PVI(i,j-1) - beta
1377	SVIMJ(I,J) =-PVI(i,j-1)
1378	!!!! OPPOSY(I,J) = ROG(1.-RO/ROW)(H(I,J-1)*2)dx/2.
1379	OPPOSY(i,j)= ROG(H(I,J-1)2)dx/2. - ROWG((min(slv(i,j-1)-B(i,j-1),0.))2)dx/2.
1380	else
1381	beta=FROTMY(I,J)*dx
1382	beta=min(beta,0.3*PVI(i,j))
1383	! SVIJ(i,j) =-6*PVI(i,j) + beta
1384	! SVIPJ(i,j) = 6*PVI(i,j)
1385	SVIJ(i,j) =-PVI(i,j) + beta
1386	SVIPJ(i,j) = PVI(i,j)
1387	!!!! OPPOSY(i,j)= ROG(1.-RO/ROW)(H(I,J)*2)dx/2.
1388	OPPOSY(i,j)= ROG(H(I,J)2)dx/2. - ROWG((min(slv(i,j)-B(i,j),0.))2)dx/2.
1389	endif
1390	END SUBROUTINE V_COIN
1391	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1392	SUBROUTINE U_CISAILL(DIR)
1393	!Rempli la matice pour le calcul de Uij
1394	IMPLICIT NONE
1395
1396	INTEGER DIR!Direction du front vers l'aval+1,vers l'amont-1
1397	! Equivaut frontfacey
1398
1399	IF (DIR==1) THEN
1400	if (.not.(flotmy(i,j).and.flotmy(i-1,j))) then
1401	! if (.not.(flgzmy(i,j).and.flgzmy(i-1,j))) then
1402	!On ne calcule pas le cisaillement sur le bord des streams
1403	if (ice(i+1,j)==0.and.ice(i-2,j)==0) then
1404	TUIJ(I,J) = 1
1405	TUIMJ(I,J) =-1
1406	else
1407	TUIJ(I,J) = 1
1408	TUPIJ(I,J) =-0.5
1409	TUMIJ(I,J) =+0.5
1410	endif
1411	else
1412	!print*, 'U_CISAILL(+1)',i,j
1413	!Le front est en haut, on calcule EXY=0 vers le bas (en arriere)
1414	!du/dy
1415	TUIJ(I,J) = 1
1416	TUIMJ(I,J) =-1
1417	!dv/dx
1418	TVIPJ(I,J) = 1
1419	TVMIPJ(I,J)=-1
1420	endif
1421	ELSE
1422	if (.not.(flotmy(i,j+1).and.flotmy(i-1,j+1))) then
1423	! if (.not.(flgzmy(i,j+1).and.flgzmy(i-1,j+1))) then
1424	!On ne calcule pas le cisaillement sur le bord des streams
1425	if (ice(i+1,j)==0.and.ice(i-2,j)==0) then
1426	TUIJ(I,J) = 1
1427	TUIPJ(I,J) =-1
1428	else
1429	TUIJ(I,J) = 1
1430	TUPIJ(I,J) =-0.5
1431	TUMIJ(I,J) =+0.5
1432	endif
1433	else
1434	!Le front est en bas, on calcule EXY=0 vers le haut(en arriere)
1435	!du/dy
1436	TUIJ(I,J) =-1
1437	TUIPJ(I,J) = 1
1438	!dv/dx
1439	TVIJ(I,J) = 1
1440	TVMIJ(I,J) =-1
1441	endif
1442	ENDIF
1443	OPPOSX(I,J)= 0
1444	END SUBROUTINE U_CISAILL
1445	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1446	SUBROUTINE V_CISAILL(DIR)
1447	!Rempli la matice pour le calcul de Vij
1448	IMPLICIT NONE
1449
1450	INTEGER DIR!Direction du front vers l'aval+1,vers l'amont-1
1451	! Equivaut frontfacex
1452
1453	IF (DIR==1) THEN
1454	! if (i.eq.54) print*,i,j,flotmx(i,j),flotmx(i,j-1)
1455	if ((.not.flotmy(i,j)).or..not.(flotmx(i,j).and.flotmx(i,j-1))) then
1456	! if (.not.(flgzmx(i,j).and.flgzmx(i,j-1))) then
1457	!On ne calcule pas le cisaillement sur le bord des streams
1458	if (ice(i,j+1)==0.and.ice(i,j-2)==0) then
1459	SVIJ(I,J) = 1
1460	SVMIJ(I,J) =-1
1461	else
1462	SVIJ(I,J) = 1
1463	SVIPJ(I,J) =-0.5
1464	SVIMJ(I,J) =+0.5
1465	endif
1466	else
1467	!Le front est a droite, on calcule EXY=0 vers la gauche (en arriere)
1468	!dv/dx
1469	SVIJ(I,J) = 1
1470	SVMIJ(I,J) =-1
1471	!du/dy
1472	SUPIJ(I,J) = 1
1473	SUPIMJ(I,J) =-1
1474	endif
1475	ELSE
1476	if ((.not.flotmy(i,j)).or..not.(flotmx(i+1,j).and.flotmx(i+1,j-1))) then
1477	! if (.not.(flgzmx(i+1,j).and.flgzmx(i+1,j-1))) then
1478	!On ne calcule pas le cisaillement sur le bord des streams
1479	if (ice(i,j+1)==0.and.ice(i,j-2)==0) then
1480	SVIJ(I,J) = 1
1481	SVPIJ(I,J) =-1
1482	else
1483	SVIJ(I,J) = 1
1484	SVIPJ(I,J) =-0.5
1485	SVIMJ(I,J) =+0.5
1486	endif
1487	else
1488	!Le front est a gauche, on calcule EXY=0 vers la droite (en arriere)
1489	!dv/dx
1490	SVIJ(I,J) =-1
1491	SVPIJ(I,J) = 1
1492	!du/dy
1493	SUIJ(I,J) = 1
1494	SUIMJ(I,J) =-1
1495	endif
1496	ENDIF
1497	OPPOSY(I,J)= 0
1498	END SUBROUTINE V_CISAILL
1499	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1500
1501	SUBROUTINE V_TONGUE
1502
1503	!la zone a 2 voisins alignes selon l'axe Y : on traite Vij
1504	beta=FROTMY(I,J)dxdx
1505	SVIJ(I,J) =-6PVI(i,j-1) - 6PVI(i,j) - beta
1506	SVIMJ(I,J) = 6*PVI(i,j-1)
1507	SVIPJ(I,J) = 6*PVI(i,j)
1508	moteur= ROGHMX(I,J) &
1509	(S(I,J)-S(I,J-1))/dx
1510	moteur=min(moteur,moteurmax)
1511	moteur=max(moteur,-moteurmax)
1512	OPPOSY(I,J)=moteurdxdx
1513	END SUBROUTINE V_TONGUE
1514
1515	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1516
1517	SUBROUTINE V_TONGUEFRONT(cas)
1518	!la langue de glace est verticale (selon les Y): calcul Vij sur le front
1519	integer CAS !valeur de frontfacey : si le front est au-dessus ou
1520	!au-dessous du point en glace
1521	if (cas==1) then!le point en glace est dessous: ICE(i,j-1)=1
1522
1523	beta=FROTMY(I,J)*dx
1524	SVIJ(I,J) = 6*PVI(i,j-1) - beta
1525	SVIMJ(I,J) =-6*PVI(i,j-1)
1526	! SVIJ(I,J) = PVI(i,j-1) - beta
1527	! SVIMJ(I,J) =-PVI(i,j-1)
1528	!!!! OPPOSY(I,J) = ROG(1.-RO/ROW)(H(I,J-1)*2)dx/2.
1529	OPPOSY(i,j)= ROG(H(I,J-1)2)dx/2. - ROWG((min(slv(i,j-1)-B(i,j-1),0.))2)dx/2.
1530
1531	else !cas=-1 le point en glace est dessus : ICE(i,j)=1
1532
1533	beta=FROTMY(I,J)*dx
1534	SVIJ(i,j) =-6*PVI(i,j) - beta
1535	SVIPJ(i,j) = 6*PVI(i,j)
1536	! SVIJ(i,j) =-PVI(i,j) - beta
1537	! SVIPJ(i,j) = PVI(i,j)
1538	!!!! OPPOSY(i,j)= ROG(1.-RO/ROW)(H(I,J)*2)dx/2.
1539	OPPOSY(i,j)= ROG(H(I,J)2)dx/2.- ROWG((min(slv(i,j)-B(i,j),0.))2)dx/2.
1540
1541	endif
1542	END SUBROUTINE V_TONGUEFRONT
1543
1544	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1545	SUBROUTINE U_TONGUE
1546
1547	!la zone a 2 voisins alignes selon l'axe X : on traite Uij
1548	beta=FROTMX(i,j)dxdx
1549	TUIJ(I,J) =-6PVI(i-1,j) - 6PVI(i,j) - beta
1550	TUMIJ(I,J) = 6*PVI(i-1,j)
1551	TUPIJ(I,J) = 6*PVI(i,j)
1552	moteur= ROGHMX(I,J)* &
1553	(S(I,J)-S(I-1,J))/dx
1554	moteur=min(moteur,moteurmax)
1555	moteur=max(moteur,-moteurmax)
1556	OPPOSX(I,J)=moteurdxdx
1557	END SUBROUTINE U_TONGUE
1558
1559	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1560
1561	SUBROUTINE U_TONGUEFRONT(cas)
1562	!langue de glace horizontale (selon les X): calcul Uij sur le front
1563	integer CAS !valeur de frontfacex :si le front est a-droite ou
1564	!a-gauche du point en glace
1565	if (cas==1) then!le point en glace est a gauche: ICE(i-1,j)=1
1566
1567	beta=FROTMX(i,j)*dx
1568	TUIJ(I,J) = 6*PVI(i-1,j) - beta
1569	TUMIJ(I,J) =-6*PVI(i-1,j)
1570	! TUIJ(I,J) = PVI(i-1,j) - beta
1571	! TUMIJ(I,J) =-PVI(i-1,j)
1572	!!!! OPPOSX(I,J) = ROG(1.-RO/ROW)(H(I-1,J)*2)dx/2.
1573	OPPOSX(I,J) = ROG(H(I-1,J)2)dx/2. - ROWG((min(slv(i-1,j)-B(i-1,j),0.))2)dx/2.
1574
1575	else !cas=-1 le point en glace est a droite: ICE(i,j)=1
1576
1577	beta=FROTMX(i,j)*dx
1578	TUIJ(I,J) =-6*PVI(i,j) - beta
1579	TUPIJ(I,J) = 6*PVI(i,j)
1580	! TUIJ(I,J) =-PVI(i,j) - beta
1581	! TUPIJ(I,J) = PVI(i,j)
1582	!!!! OPPOSX(I,J) = ROG(1.-RO/ROW)(H(I,J)*2)dx/2.
1583	OPPOSX(I,J) = ROG(H(I,J)2)dx/2. - ROWG((min(slv(i,j)-B(i,j),0.))2)dx/2.
1584
1585	endif
1586	END SUBROUTINE U_TONGUEFRONT
1587
1588	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1589
1590	SUBROUTINE V_TONGUE_UP(front)
1591	integer front ! valeur de frontfacey(i,j-1)
1592	OPPOSY(I,J)= 0
1593	if (front==0) then!voisin en glace (en bas) a gauche et a doite
1594	SVIJ(I,J) = 1
1595	SVMIJ(I,J) =-.25
1596	SVMIMJ(I,J)=-.25
1597	SVPIJ(I,J) =-.25
1598	SVPIMJ(I,J)=-.25
1599	SUPIMJ(I,J)=-.5
1600	SUIMJ(I,J) = .5
1601	elseif (front==1) then!voisin en glace (en bas) a gauche
1602	SVIJ(I,J) = 1
1603	SVMIJ(I,J) =-.5
1604	SVMIMJ(I,J)=-.5
1605	SUPIMJ(I,J)=-.5
1606	SUIMJ(I,J) = .5
1607	else !!voisin en glace (en bas) a doite
1608	SVIJ(I,J) = 1
1609	SVPIJ(I,J) =-.5
1610	SVPIMJ(I,J)=-.5
1611	SUPIMJ(I,J)=-.5
1612	SUIMJ(I,J) = .5
1613	endif
1614	END SUBROUTINE V_TONGUE_UP
1615
1616	SUBROUTINE V_TONGUE_DO(front)!DO=down
1617	integer front ! valeur de frontfacey(i,j)
1618
1619	OPPOSY(I,J)= 0
1620	if (front==0) then!voisin en glace a gauche et a doite
1621	SVIJ(I,J) = 1
1622	SVMIPJ(I,J)=-.25
1623	SVMIJ(I,J) =-.25
1624	SVPIPJ(I,J)=-.25
1625	SVPIJ(I,J) =-.25
1626	SUPIJ(I,J) = .5
1627	SUIJ(I,J) =-.5
1628	elseif (front==1) then!voisin en glace a gauche
1629	SVIJ(I,J) = 1
1630	SVMIPJ(I,J)=-.5
1631	SVMIJ(I,J) =-.5
1632	SUPIJ(I,J) = .5
1633	SUIJ(I,J) =-.5
1634	else!voisin en glace a doite
1635	SVIJ(I,J) = 1
1636	SVPIPJ(I,J)=-.5
1637	SVPIJ(I,J) =-.5
1638	SUPIJ(I,J) = .5
1639	SUIJ(I,J) =-.5
1640	endif
1641	END SUBROUTINE V_TONGUE_DO
1642
1643	SUBROUTINE U_TONGUE_LE(front)!LE=left
1644	integer front ! valeur de frontfacex(i,j)
1645
1646	OPPOSX(I,J)= 0
1647	if (front==0) then!voisin en glace en haut et en bas
1648	TUIJ(I,J) = 1
1649	TUPIMJ(I,J)=-.25
1650	TUIMJ(I,J) =-.25
1651	TUPIPJ(i,j)=-.25
1652	TUIPJ(i,j) =-.25
1653	TVIPJ(I,J) = .5
1654	TVIJ(I,J) =-.5
1655	elseif(front==1) then!voisin en glace en bas
1656	TUIJ(I,J) = 1
1657	TUPIMJ(I,J)=-.5
1658	TUIMJ(I,J) =-.5
1659	TVIPJ(I,J) = .5
1660	TVIJ(I,J) =-.5
1661	else!voisin en glace en haut
1662	TUIJ(i,j) =1
1663	TUIPJ(i,j) =-0.5
1664	TUPIPJ(i,j)=-0.5
1665	TVIPJ(i,j) = 0.5
1666	TVIJ(i,j) =-0.5
1667	endif
1668	END SUBROUTINE U_TONGUE_LE
1669
1670	SUBROUTINE U_TONGUE_RI(front)
1671	integer front ! valeur de frontfacex(i-1,j)
1672
1673	OPPOSX(I,J)= 0
1674	if (front==0) then!voisin en glace (a gauche) en haut et en bas
1675	TUIJ(I,J) = 1
1676	TUIMJ(I,J) =-.25
1677	TUMIMJ(I,J)=-.25
1678	TUMIPJ(i,j)=-.25
1679	TUIPJ(i,j) =-.25
1680	TVMIPJ(I,J)=-.5
1681	TVMIJ(I,J) = .5
1682	elseif (front==1) then!voisin en glace (a gauche) en bas
1683	TUIJ(I,J) = 1
1684	TUIMJ(I,J) =-.5
1685	TUMIMJ(I,J)=-.5
1686	TVMIPJ(I,J)=-.5
1687	TVMIJ(I,J) = .5
1688	else!!voisin en glace (a gauche) en haut
1689	TUIJ(i,j) =1
1690	TUMIPJ(i,j)=-0.5
1691	TUIPJ(i,j) =-0.5
1692	TVMIPJ(i,j)=-0.5
1693	TVMIJ(i,j) = 0.5
1694	endif
1695	END SUBROUTINE U_TONGUE_RI
1696	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1697	! --------------------------------------------------------------
1698	! subroutine de calcul de la largeur de bande
1699	! version sous-domaine
1700	! --------------------------------------------------------------
1701
1702	subroutine larg_bande(ksur,ksous)
1703
1704	implicit none
1705	INTEGER KSUR,KSOUS,KMIN,KMAX,IKSUR,IKSOUS,JKSUR,JKSOUS,ktest
1706	KSOUS=-1
1707	KSUR=-1
1708
1709	! noeuds U
1710
1711	! do I=1,NX
1712	! do J=1,NY
1713	do I=2,NX-1
1714	do J=2,NY-1
1715
1716	if (OKUMAT(I,J)) then
1717
1718	kmin=ligu(i,j)
1719	kmax=ligu(i,j)
1720
1721	ktest=ligu(i-1,j-1) !Ui-1
1722	if (ktest.gt.0) then
1723	kmin=min(ktest,kmin)
1724	kmax=max(ktest,kmax)
1725	endif
1726
1727	ktest=ligu(i-1,j)
1728	if (ktest.gt.0) then
1729	kmin=min(ktest,kmin)
1730	kmax=max(ktest,kmax)
1731	endif
1732
1733	ktest=ligu(i-1,j+1)
1734	if (ktest.gt.0) then
1735	kmin=min(ktest,kmin)
1736	kmax=max(ktest,kmax)
1737	endif
1738
1739	ktest=ligu(i,j-1) !Ui
1740	if (ktest.gt.0) then
1741	kmin=min(ktest,kmin)
1742	kmax=max(ktest,kmax)
1743	endif
1744
1745	ktest=ligu(i,j+1)
1746	if (ktest.gt.0) then
1747	kmin=min(ktest,kmin)
1748	kmax=max(ktest,kmax)
1749	endif
1750
1751	ktest=ligu(i+1,j-1) !Ui-1
1752	if (ktest.gt.0) then
1753	kmin=min(ktest,kmin)
1754	kmax=max(ktest,kmax)
1755	endif
1756
1757	ktest=ligu(i+1,j)
1758	if (ktest.gt.0) then
1759	kmin=min(ktest,kmin)
1760	kmax=max(ktest,kmax)
1761	endif
1762
1763	ktest=ligu(i+1,j+1)
1764	if (ktest.gt.0) then
1765	kmin=min(ktest,kmin)
1766	kmax=max(ktest,kmax)
1767	endif
1768
1769	ktest=ligv(i-1,j) !Vi-1
1770	if (ktest.gt.0) then
1771	kmin=min(ktest,kmin)
1772	kmax=max(ktest,kmax)
1773	endif
1774
1775	ktest=ligv(i-1,j+1)
1776	if (ktest.gt.0) then
1777	kmin=min(ktest,kmin)
1778	kmax=max(ktest,kmax)
1779	endif
1780
1781	ktest=ligv(i,j) !Vi
1782	if (ktest.gt.0) then
1783	kmin=min(ktest,kmin)
1784	kmax=max(ktest,kmax)
1785	endif
1786
1787	ktest=ligv(i,j+1)
1788	if (ktest.gt.0) then
1789	kmin=min(ktest,kmin)
1790	kmax=max(ktest,kmax)
1791	endif
1792	! ------------------------------- fin de recherche min max de la ligne
1793
1794	if (LIGU(I,J)-KMIN.gt.KSOUS) then
1795	KSOUS=LIGU(I,J)-KMIN
1796	IKSOUS=I
1797	JKSOUS=J
1798	endif
1799
1800	if (KMAX-LIGU(I,J).gt.KSUR) then
1801	KSUR=KMAX-LIGU(I,J)
1802	IKSUR=I
1803	JKSUR=J
1804	endif
1805
1806	endif
1807
1808
1809
1810	end do
1811	end do
1812
1813	! ----------------- impression de largeur de bande -------------------
1814	! write(6,*)'largeur de bande pour les noeuds internes en U'
1815	! write(6,*)'ksous=',ksous,' pour i=',iksous,' pour j=',jksous
1816	! write(6,*)'ksur=',ksur,' pour i=',iksur,' pour j=',jksur
1817	! --------------------------------------------------------------------
1818
1819	! noeuds V
1820	! KSOUS=-1
1821	! KSUR=-1
1822
1823	! do I=1,NX
1824	! do J=1,NY
1825	do I=2,NX-1
1826	do J=2,NY-1
1827	if (OKVMAT(I,J)) then
1828
1829	kmin=ligv(i,j)
1830	kmax=ligv(i,j)
1831
1832	ktest=ligv(i-1,j-1) !Vi-1
1833	if (ktest.gt.0) then
1834	kmin=min(ktest,kmin)
1835	kmax=max(ktest,kmax)
1836	endif
1837
1838	ktest=ligv(i-1,j)
1839	if (ktest.gt.0) then
1840	kmin=min(ktest,kmin)
1841	kmax=max(ktest,kmax)
1842	endif
1843
1844	ktest=ligv(i-1,j+1)
1845	if (ktest.gt.0) then
1846	kmin=min(ktest,kmin)
1847	kmax=max(ktest,kmax)
1848	endif
1849
1850	ktest=ligv(i,j-1) !Vi
1851	if (ktest.gt.0) then
1852	kmin=min(ktest,kmin)
1853	kmax=max(ktest,kmax)
1854	endif
1855
1856	ktest=ligv(i,j+1)
1857	if (ktest.gt.0) then
1858	kmin=min(ktest,kmin)
1859	kmax=max(ktest,kmax)
1860	endif
1861
1862	ktest=ligv(i+1,j-1) !Vi-1
1863	if (ktest.gt.0) then
1864	kmin=min(ktest,kmin)
1865	kmax=max(ktest,kmax)
1866	endif
1867
1868	ktest=ligv(i+1,j)
1869	if (ktest.gt.0) then
1870	kmin=min(ktest,kmin)
1871	kmax=max(ktest,kmax)
1872	endif
1873
1874	ktest=ligv(i+1,j+1)
1875	if (ktest.gt.0) then
1876	kmin=min(ktest,kmin)
1877	kmax=max(ktest,kmax)
1878	endif
1879
1880	ktest=ligu(i,j-1) !Ui
1881	if (ktest.gt.0) then
1882	kmin=min(ktest,kmin)
1883	kmax=max(ktest,kmax)
1884	endif
1885
1886	ktest=ligu(i,j)
1887	if (ktest.gt.0) then
1888	kmin=min(ktest,kmin)
1889	kmax=max(ktest,kmax)
1890	endif
1891
1892	ktest=ligu(i+1,j-1) !Ui-1
1893	if (ktest.gt.0) then
1894	kmin=min(ktest,kmin)
1895	kmax=max(ktest,kmax)
1896	endif
1897
1898	ktest=ligu(i+1,j)
1899	if (ktest.gt.0) then
1900	kmin=min(ktest,kmin)
1901	kmax=max(ktest,kmax)
1902	endif
1903	! ------------------------------- fin de recherche min max de la ligne
1904	if (LIGV(I,J)-KMIN.gt.KSOUS) then
1905	KSOUS=LIGV(I,J)-KMIN
1906	IKSOUS=I
1907	JKSOUS=J
1908	endif
1909
1910	if (KMAX-LIGV(I,J).gt.KSUR) then
1911	KSUR=KMAX-LIGV(I,J)
1912	IKSUR=I
1913	JKSUR=J
1914	endif
1915
1916	endif
1917	end do
1918	end do
1919
1920	! ----------------- impression de largeur de bande -------------------
1921	! write(6,*)'largeur de bande pour les noeuds internes en V'
1922	! write(6,*)'ksous=',ksous,' pour i=',iksous,' pour j=',jksous
1923	! write(6,*)'ksur=',ksur,' pour i=',iksur,' pour j=',jksur
1924	! ---------------------------------------------------------------
1925
1926	end subroutine larg_bande
1927	!----------------------------------------------------------------------
1928	! FIN DU CALCUL DE LARGEUR DE BANDE
1929	!----------------------------------------------------------------------
1930
1931	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1932
1933	SUBROUTINE DEBUGAGE
1934	implicit none
1935	logical yesno
1936	! character aaa
1937	integer iii,jjj
1938	character(len=80) :: filin
1939	real,dimension(-1:1,-1:1) :: Au, Bv
1940	real,dimension(-1:0,0:1) :: Av
1941	real,dimension(0:1,-1:0) :: Bu
1942	yesno=.false.
1943	if (yesno) then
1944	do i=2,nx-1
1945	do j=2,ny-1
1946	Au(-1,-1)=TUMIMJ(I,J)
1947	Au(-1,0)=TUMIJ(I,J)
1948	Au(-1,1)=TUMIPJ(I,J)
1949	Au(0,-1)=TUIMJ(I,J)
1950	Au(0,0)=TUIJ(I,J)
1951	Au(0,1)=TUIPJ(I,J)
1952	Au(1,-1)=TUPIMJ(I,J)
1953	Au(1,0)=TUPIJ(I,J)
1954	Au(1,1)=TUPIPJ(I,J)
1955	!!!!!!!!
1956	Av(-1,0)=TVMIJ(I,J)
1957	Av(-1,1)=TVMIPJ(I,J)
1958	Av(0,0)=TVIJ(I,J)
1959	Av(0,1)=TVIPJ(I,J)
1960	!!!!!!!!
1961	Bv(-1,-1)=SVMIMJ(I,J)
1962	Bv(-1,0)=SVMIJ(I,J)
1963	Bv(-1,1)=SVMIPJ(I,J)
1964	Bv(0,-1)=SVIMJ(I,J)
1965	Bv(0,0)=SVIJ(I,J)
1966	Bv(0,1)=SVIPJ(I,J)
1967	Bv(1,-1)=SVPIMJ(I,J)
1968	Bv(1,0)=SVPIJ(I,J)
1969	Bv(1,1)=SVPIPJ(I,J)
1970	!!!!!!!!
1971	Bu(0,0)=SUIJ(I,J)
1972	Bu(0,-1)=SUIMJ(I,J)
1973	Bu(1,0)=SUPIJ(I,J)
1974	Bu(1,-1)=SUPIMJ(I,J)
1975	do iii=-1,1
1976	do jjj=-1,1
1977	if ((Au(iii,jjj).ne.0).and. &
1978	(abs(OPPOSX(I+iii,J+jjj)/Au(0,0)).gt.100)) then
1979	write(6,*) 'erreur i,j=',i,j,'Au(',iii,jjj,')'
1980	write(6,*) 'OPPOSX=',OPPOSX(I+iii,J+jjj),'methode',METHOD(I,J,1),METHOD(I,J,2)
1981	! if(iii.ne.0.or.jjj.ne.0)read*, aaa!stop
1982	write(6,*)'Au(0,0)=',Au(0,0),'OPPOSX/Au(0,0)',OPPOSX(I+iii,J+jjj)/Au(0,0)
1983	endif
1984	if ((Bv(iii,jjj).ne.0).and. &
1985	(abs(OPPOSY(I+iii,J+jjj)/Bv(0,0)).gt.100)) then
1986	write(6,*) 'erreuri,j=',i,j,'Bv(',iii,jjj,')'
1987	write(6,*) 'OPPOSY=',OPPOSY(I+iii,J+jjj),'methode',METHOD(I,J,1),METHOD(I,J,2)
1988	write(6,*)'Bv(0,0)=',Bv(0,0),'OPPOSY/Bv(0,0)',OPPOSY(I+iii,J+jjj)/Bv(0,0)
1989	! if(iii.ne.0.or.jjj.ne.0)read*, aaa!stop
1990	endif
1991	enddo
1992	enddo
1993	do iii=-1,0
1994	do jjj=0,1
1995	if ((Av(iii,jjj).ne.0).and. &
1996	(abs(OPPOSY(I+iii,J+jjj)/Au(0,0)).gt.100)) then
1997	write(6,*) 'erreur i,j=',i,j,'Av(',iii,jjj,')'
1998	write(6,*) 'OPPOSY=',OPPOSY(I+iii,J+jjj),'methode',METHOD(I,J,1),METHOD(I,J,2)
1999	write(6,*) 'AV(0,0)=',Av(0,0),'OPPOSY/Av(0,0)',OPPOSY(I+iii,J+jjj)/Av(0,0)
2000
2001	! if(iii.ne.0.or.jjj.ne.0)read*, aaa!stop
2002	endif
2003	enddo
2004	enddo
2005	do jjj=-1,0
2006	do iii=0,1
2007	if ((Bu(iii,jjj).ne.0).and. &
2008	(abs(OPPOSX(I+iii,J+jjj)/Bv(0,0)).gt.100)) then
2009	write(6,*) 'erreur i,j=',i,j,'Bu(',iii,jjj,')'
2010	write(6,*) 'OPPOSX=',OPPOSX(I+iii,J+jjj),'methode',METHOD(I,J,1),METHOD(I,J,2)
2011	write(6,*) 'Bu(0,0)=',Bu(0,0),'OPPOSX/Bu(0,0)',OPPOSX(I+iii,J+jjj)/Bu(0,0)
2012	! if(iii.ne.0.or.jjj.ne.0) read*, aaa!stop
2013	endif
2014	enddo
2015	enddo
2016	write(6,*)'==================================================='
2017	enddo
2018	enddo
2019	endif
2020	!!!!!!!!!!!impression de la methode!!!!!!!!!!!!!!!!!
2021	filin = 'METHOD_uv'
2022
2023	open (UNIT=num_file3,file=filin)
2024	write(num_file3,*) geoplace,'time=',time
2025	write(num_file3,) NXNY,DX,NX,NY
2026	do j=1,ny
2027	do i=1,nx
2028	! do i=35,60
2029	! do j=35,60
2030	write(num_file3,*) i,j, METHOD(i,j,1),METHOD(i,j,2)
2031	enddo
2032	enddo
2033	close (num_file3)
2034	filin = 'TABLE_met_U'
2035
2036	open (UNIT=num_file3,file=filin)
2037	write(num_file3,*) geoplace,'time=',time
2038	write(num_file3,) NXNY,DX,NX,NY
2039	do j=1,ny
2040	do i=1,nx
2041	write(num_file3,*) i,j, METHOD2(i,j,1)
2042	enddo
2043	enddo
2044	close (num_file3)
2045	filin = 'TABLE_met_V'
2046
2047	open (UNIT=num_file3,file=filin)
2048	write(num_file3,*) geoplace,'time=',time
2049	write(num_file3,) NXNY,DX,NX,NY
2050	do j=1,ny
2051	do i=1,nx
2052	write(num_file3,*) i,j, METHOD2(i,j,2)
2053	enddo
2054	enddo
2055	close (num_file3)
2056	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2057	!!!!!!!!!!!impression tuij!!!!!!!!!!!!!!!!!
2058	filin = 'TABLE_tuij'
2059
2060	open(UNIT=num_file4,file=filin)
2061	write(num_file4,*) geoplace,'time=',time
2062	write(num_file4,) NXNY,DX,NX,NY
2063	do i=1,nx
2064	do j=1,ny
2065	write(num_file4,*) i,j,tuij(i,j)!/tuij(i,j)
2066	enddo
2067	enddo
2068	close (num_file4)
2069	filin = 'TABLE_tumimj'
2070
2071	open(UNIT=num_file4,file=filin)
2072	write(num_file4,*) geoplace,'time=',time
2073	write(num_file4,) NXNY,DX,NX,NY
2074	do i=1,nx
2075	do j=1,ny
2076	write(num_file4,*) i,j,tumimj(i,j)/tuij(i,j)
2077	enddo
2078	enddo
2079	close (num_file4)
2080	filin = 'TABLE_tumij'
2081
2082	open(UNIT=num_file4,file=filin)
2083	write(num_file4,*) geoplace,'time=',time
2084	write(num_file4,) NXNY,DX,NX,NY
2085	do i=1,nx
2086	do j=1,ny
2087	write(num_file4,*) i,j,tumij(i,j)/tuij(i,j)
2088	enddo
2089	enddo
2090	close (num_file4)
2091	filin = 'TABLE_tuimj'
2092
2093	open(UNIT=num_file4,file=filin)
2094	write(num_file4,*) geoplace,'time=',time
2095	write(num_file4,) NXNY,DX,NX,NY
2096	do i=1,nx
2097	do j=1,ny
2098	write(num_file4,*) i,j,tuimj(i,j)/tuij(i,j)
2099	enddo
2100	enddo
2101	close (num_file4)
2102	filin = 'TABLE_tumipj'
2103
2104	open(UNIT=num_file4,file=filin)
2105	write(num_file4,*) geoplace,'time=',time
2106	write(num_file4,) NXNY,DX,NX,NY
2107	do i=1,nx
2108	do j=1,ny
2109	write(num_file4,*) i,j,tumipj(i,j)/tuij(i,j)
2110	enddo
2111	enddo
2112	close (num_file4)
2113	filin = 'TABLE_tuimj'
2114
2115	open(UNIT=num_file4,file=filin)
2116	write(num_file4,*) geoplace,'time=',time
2117	write(num_file4,) NXNY,DX,NX,NY
2118	do i=1,nx
2119	do j=1,ny
2120	write(num_file4,*) i,j,tuimj(i,j)/tuij(i,j)
2121	enddo
2122	enddo
2123	close (num_file4)
2124	filin = 'TABLE_tuipj'
2125
2126	open(UNIT=num_file4,file=filin)
2127	write(num_file4,*) geoplace,'time=',time
2128	write(num_file4,) NXNY,DX,NX,NY
2129	do i=1,nx
2130	do j=1,ny
2131	write(num_file4,*) i,j,tuipj(i,j)/tuij(i,j)
2132	enddo
2133	enddo
2134	close (num_file4)
2135	filin = 'TABLE_tupimj'
2136
2137	open(UNIT=num_file4,file=filin)
2138	write(num_file4,*) geoplace,'time=',time
2139	write(num_file4,) NXNY,DX,NX,NY
2140	do i=1,nx
2141	do j=1,ny
2142	write(num_file4,*) i,j,tupimj(i,j)/tuij(i,j)
2143	enddo
2144	enddo
2145	close (num_file4)
2146	filin = 'TABLE_tupij'
2147
2148	open(UNIT=num_file4,file=filin)
2149	write(num_file4,*) geoplace,'time=',time
2150	write(num_file4,) NXNY,DX,NX,NY
2151	do i=1,nx
2152	do j=1,ny
2153	write(num_file4,*) i,j,tupij(i,j)/tuij(i,j)
2154	enddo
2155	enddo
2156	close (num_file4)
2157	filin = 'TABLE_tupipj'
2158
2159	open(UNIT=num_file4,file=filin)
2160	write(num_file4,*) geoplace,'time=',time
2161	write(num_file4,) NXNY,DX,NX,NY
2162	do i=1,nx
2163	do j=1,ny
2164	write(num_file4,*) i,j,tupipj(i,j)/tuij(i,j)
2165	enddo
2166	enddo
2167	close (num_file4)
2168	filin = 'TABLE_tvij'
2169
2170	open(UNIT=num_file4,file=filin)
2171	write(num_file4,*) geoplace,'time=',time
2172	write(num_file4,) NXNY,DX,NX,NY
2173	do i=1,nx
2174	do j=1,ny
2175	write(num_file4,*) i,j,tvij(i,j)/tuij(i,j)
2176	enddo
2177	enddo
2178	close (num_file4)
2179	filin = 'TABLE_tvmij'
2180
2181	open(UNIT=num_file4,file=filin)
2182	write(num_file4,*) geoplace,'time=',time
2183	write(num_file4,) NXNY,DX,NX,NY
2184	do i=1,nx
2185	do j=1,ny
2186	write(num_file4,*) i,j,tvmij(i,j)/tuij(i,j)
2187	enddo
2188	enddo
2189	close (num_file4)
2190	filin = 'TABLE_tvmipj'
2191
2192	open(UNIT=num_file4,file=filin)
2193	write(num_file4,*) geoplace,'time=',time
2194	write(num_file4,) NXNY,DX,NX,NY
2195	do i=1,nx
2196	do j=1,ny
2197	write(num_file4,*) i,j,tvmipj(i,j)/tuij(i,j)
2198	enddo
2199	enddo
2200	close (num_file4)
2201	filin = 'TABLE_tvipj'
2202
2203	open(UNIT=num_file4,file=filin)
2204	write(num_file4,*) geoplace,'time=',time
2205	write(num_file4,) NXNY,DX,NX,NY
2206	do i=1,nx
2207	do j=1,ny
2208	write(num_file4,*) i,j,tvipj(i,j)/tuij(i,j)
2209	enddo
2210	enddo
2211	close (num_file4)
2212	filin = 'TABLE_opposx'
2213
2214	open(UNIT=num_file4,file=filin)
2215	write(num_file4,*) geoplace,'time=',time
2216	write(num_file4,) NXNY,DX,NX,NY
2217	do j=1,ny
2218	do i=1,nx
2219	write(num_file4,*) i,j,opposx(i,j)/tuij(i,j)
2220	enddo
2221	enddo
2222	close (num_file4)
2223	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2224	!!!!!!!!!!!impression svij!!!!!!!!!!!!!!!!!
2225	filin = 'TABLE_svij'
2226
2227	open(UNIT=num_file4,file=filin)
2228	write(num_file4,*) geoplace,'time=',time
2229	write(num_file4,) NXNY,DX,NX,NY
2230	do i=1,nx
2231	do j=1,ny
2232	write(num_file4,*) i,j,svij(i,j)!/svij(i,j)
2233	enddo
2234	enddo
2235	close (num_file4)
2236	filin = 'TABLE_svmimj'
2237
2238	open(UNIT=num_file4,file=filin)
2239	write(num_file4,*) geoplace,'time=',time
2240	write(num_file4,) NXNY,DX,NX,NY
2241	do i=1,nx
2242	do j=1,ny
2243	write(num_file4,*) i,j,svmimj(i,j)/svij(i,j)
2244	enddo
2245	enddo
2246	close (num_file4)
2247	filin = 'TABLE_svimj'
2248
2249	open(UNIT=num_file4,file=filin)
2250	write(num_file4,*) geoplace,'time=',time
2251	write(num_file4,) NXNY,DX,NX,NY
2252	do i=1,nx
2253	do j=1,ny
2254	write(num_file4,*) i,j,svimj(i,j)/svij(i,j)
2255	enddo
2256	enddo
2257	close (num_file4)
2258	filin = 'TABLE_svmipj'
2259
2260	open(UNIT=num_file4,file=filin)
2261	write(num_file4,*) geoplace,'time=',time
2262	write(num_file4,) NXNY,DX,NX,NY
2263	do i=1,nx
2264	do j=1,ny
2265	write(num_file4,*) i,j,svmipj(i,j)/svij(i,j)
2266	enddo
2267	enddo
2268	close (num_file4)
2269	filin = 'TABLE_svimj'
2270
2271	open(UNIT=num_file4,file=filin)
2272	write(num_file4,*) geoplace,'time=',time
2273	write(num_file4,) NXNY,DX,NX,NY
2274	do i=1,nx
2275	do j=1,ny
2276	write(num_file4,*) i,j,svimj(i,j)/svij(i,j)
2277	enddo
2278	enddo
2279	close (num_file4)
2280	filin = 'TABLE_svipj'
2281
2282	open(UNIT=num_file4,file=filin)
2283	write(num_file4,*) geoplace,'time=',time
2284	write(num_file4,) NXNY,DX,NX,NY
2285	do i=1,nx
2286	do j=1,ny
2287	write(num_file4,*) i,j,svipj(i,j)/svij(i,j)
2288	enddo
2289	enddo
2290	close (num_file4)
2291	filin = 'TABLE_svpimj'
2292
2293	open(UNIT=num_file4,file=filin)
2294	write(num_file4,*) geoplace,'time=',time
2295	write(num_file4,) NXNY,DX,NX,NY
2296	do i=1,nx
2297	do j=1,ny
2298	write(num_file4,*) i,j,svpimj(i,j)/svij(i,j)
2299	enddo
2300	enddo
2301
2302	close (num_file4)
2303	filin = 'TABLE_svpij'
2304
2305	open(UNIT=num_file4,file=filin)
2306	write(num_file4,*) geoplace,'time=',time
2307	write(num_file4,) NXNY,DX,NX,NY
2308	do i=1,nx
2309	do j=1,ny
2310	write(num_file4,*) i,j,svpij(i,j)/svij(i,j)
2311	enddo
2312	enddo
2313	close (num_file4)
2314	filin = 'TABLE_svpipj'
2315
2316	open(UNIT=num_file4,file=filin)
2317	write(num_file4,*) geoplace,'time=',time
2318	write(num_file4,) NXNY,DX,NX,NY
2319	do i=1,nx
2320	do j=1,ny
2321	write(num_file4,*) i,j,svpipj(i,j)/svij(i,j)
2322	enddo
2323	enddo
2324	close (num_file4)
2325	filin = 'TABLE_suij'
2326
2327	open(UNIT=num_file4,file=filin)
2328	write(num_file4,*) geoplace,'time=',time
2329	write(num_file4,) NXNY,DX,NX,NY
2330	do i=1,nx
2331	do j=1,ny
2332	write(num_file4,*) i,j,suij(i,j)/svij(i,j)
2333	enddo
2334	enddo
2335	close (num_file4)
2336	filin = 'TABLE_suimj'
2337
2338	open(UNIT=num_file4,file=filin)
2339	write(num_file4,*) geoplace,'time=',time
2340	write(num_file4,) NXNY,DX,NX,NY
2341	do i=1,nx
2342	do j=1,ny
2343	write(num_file4,*) i,j,suimj(i,j)/svij(i,j)
2344	enddo
2345	enddo
2346	close (num_file4)
2347	filin = 'TABLE_supimj'
2348
2349	open(UNIT=num_file4,file=filin)
2350	write(num_file4,*) geoplace,'time=',time
2351	write(num_file4,) NXNY,DX,NX,NY
2352	do i=1,nx
2353	do j=1,ny
2354	write(num_file4,*) i,j,supimj(i,j)/svij(i,j)
2355	enddo
2356	enddo
2357	close (num_file4)
2358	filin = 'TABLE_supij'
2359
2360	open(UNIT=num_file4,file=filin)
2361	write(num_file4,*) geoplace,'time=',time
2362	write(num_file4,) NXNY,DX,NX,NY
2363	do i=1,nx
2364	do j=1,ny
2365	write(num_file4,*) i,j,supij(i,j)/svij(i,j)
2366	enddo
2367	enddo
2368	close (num_file4)
2369	filin = 'TABLE_opposy'
2370
2371	open(UNIT=num_file4,file=filin)
2372	write(num_file4,*) geoplace,'time=',time
2373	write(num_file4,) NXNY,DX,NX,NY
2374	do j=1,ny
2375	do i=1,nx
2376	write(num_file4,*) i,j,opposy(i,j)/svij(i,j)
2377	enddo
2378	enddo
2379	close (num_file4)
2380	filin = 'TABLE_FRONT'
2381
2382	open(UNIT=num_file4,file=filin)
2383	write(num_file4,*) geoplace,'time=',time
2384	write(num_file4,) NXNY,DX,NX,NY
2385	do i=1,nx
2386	do j=1,ny
2387	IF (front(I,J).EQ.4) THEN
2388	INEARX=IFLOTMX(I-1,J) +IFLOTMX(I,J) +IFLOTMX(I+1,J) &
2389	+IFLOTMX(I-1,J-1)+IFLOTMX(I,J-1)+IFLOTMX(I+1,J-1) &
2390	+IFLOTMX(I-1,J+1)+IFLOTMX(I,J+1)+IFLOTMX(I+1,J+1) &
2391	+IFLOTMY(I-1,J) +IFLOTMY(I,J) +IFLOTMY(I+1,J) &
2392	+IFLOTMY(I-1,J-1)+IFLOTMY(I,J-1)+IFLOTMY(I+1,J-1) &
2393	+IFLOTMY(I-1,J+1)+IFLOTMY(I,J+1)+IFLOTMY(I+1,J+1)
2394	if (INEARX==0) then
2395	write(num_file4,*) i,j,' 6'
2396	else
2397	write(num_file4,*) i,j,front(i,j)
2398	endif
2399	ELSE
2400	write(num_file4,*) i,j,front(i,j)
2401	INEARX=1000
2402	END IF
2403
2404	enddo
2405	enddo
2406	close (num_file4)
2407	filin = 'sgbsv_u.dat'
2408
2409	open(UNIT=num_file4,file=filin)
2410	write(num_file4,*) geoplace,'time=',time
2411	write(num_file4,) NXNY,DX,NX,NY
2412	do i=1,nx
2413	do j=1,ny
2414	write(num_file4,*) i,j,uxbar(i,j)
2415	enddo
2416	enddo
2417	close (num_file4)
2418	filin = 'sgbsv_v.dat'
2419
2420	open(UNIT=num_file4,file=filin)
2421	write(num_file4,*) geoplace,'time=',time
2422	write(num_file4,) NXNY,DX,NX,NY
2423	do i=1,nx
2424	do j=1,ny
2425	write(num_file4,*) i,j,uybar(i,j)
2426	enddo
2427	enddo
2428	close (num_file4)
2429	filin = 'TABLE_h'
2430
2431	open(UNIT=num_file4,file=filin)
2432	write(num_file4,*) geoplace,'time=',time
2433	write(num_file4,) NXNY,DX,NX,NY
2434	do i=1,nx
2435	do j=1,ny
2436	write(num_file4,*) i,j,h(i,j)
2437	enddo
2438	enddo
2439	close (num_file4)
2440	filin = 'TABLE_ICE'
2441
2442	open(UNIT=num_file4,file=filin)
2443	write(num_file4,*) geoplace,'time=',time
2444	write(num_file4,) NXNY,DX,NX,NY
2445	do i=1,nx
2446	do j=1,ny
2447	write(num_file4,*) i,j,ICE(i,j)
2448	enddo
2449	enddo
2450	close (num_file4)
2451	filin = 'TABLE_index'
2452
2453	open(UNIT=num_file4,file=filin)
2454	write(num_file4,*) geoplace,'time=',time
2455	write(num_file4,) NXNY,DX,NX,NY
2456	do i=1,nx
2457	do j=1,ny
2458	if (okumat(i,j).or.okvmat(i,j)) then
2459	write(num_file4,*) i,j,ligu(i,j),ligv(i,j)
2460	endif
2461	enddo
2462	enddo
2463	close (num_file4)
2464	filin = 'TABLE_pvi'
2465
2466	open(UNIT=num_file4,file=filin)
2467	write(num_file4,*) geoplace,'time=',time
2468	write(num_file4,) NXNY,DX,NX,NY
2469	do i=1,nx
2470	do j=1,ny
2471	write(num_file4,*) i,j,pvi(i,j)
2472	enddo
2473	enddo
2474	close (num_file4)
2475	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2476	END SUBROUTINE DEBUGAGE
2477	! RETURN
2478	END SUBROUTINE REMPLIDOM

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

Context Navigation

source: trunk/SOURCES/Old-sources/remplimat-ant-0.5-40km.f90 @ 159

Download in other formats: