Context Navigation

source: trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90 @ 74

Last change on this file since 74 was 74, checked in by dumas, 8 years ago
OpenMP parallelization in remplimat-shelves-tabTu, litho and taubed
File size: 42.7 KB

Line
1	! La matrice du systeme elliptique est nommee "L2" dans ce qui suit
2	! mais n'est pas dimensionnee car elle serait trop grande (2 nx ny x 2 nx ny)
3
4
5	! Au lieu de cela on ecrit les elements non nuls de chaque ligne
6	! Tu et Tv pour les lignes de l'equation en U
7	! Su et Sv pour les lignes de l'equation en V
8
9
10	! dans cette version les Tuij sont sous forme de tableau.
11	! le conditions aux limites sont donnees sur les bords de la grille
12	! Il y a un nettoyage de la matrice pour enlever les noeuds qui ne servent pas.
13
14
15	!--------------------------------------------------------------------------
16	subroutine rempli_L2(nx1,nx2,ny1,ny2,uxprec,uyprec,uxnew,uynew,imx,imy,ifail_L2)
17	!
18	! - ecrit l'equation elliptique des vitesses sur le domaine nx1,nx2 et ny1,ny2
19	! - appelle la routine de resolution : call resol_ellipt
20	! - renvoie les nouvelles vitesses uxnew, uynew
21	!
22	!
23	! nx1,nx2 bornes du domaine
24	! ny1,ny2
25	!
26	! uxprex(nx1:nx2,ny1:ny2)
27	! uyprec(nx1:nx2,ny1:ny2) vitesses de l'iteration precedente
28	!
29	! uxnew(nx1:nx2,ny1:ny2)
30	! uynew(nx1:nx2,ny1:ny2) uynew resultat de cette iteration
31	!
32	! imx(nx1:nx2,ny1:ny2) masque pour imposer les vitesses ou leur dérivee
33	! imy(nx1:nx2,ny1:ny2) masque pour imposer les vitesses ou leur dérivée
34	!
35	! eventuellement le domaine n1,n2 peut etre un sous-domaine de nx,ny
36	! attention, dans ce cas l'appel devra être
37	! call rempli-L2 (nx1,nx2,ny1,ny2,ux1(nx1:nx2,ny1:ny2),uy1(nx1:nx2,ny1:ny2), &
38	! ux2(nx1:nx2,ny1:ny2),uy2(nx1:nx2,ny1:ny2), &
39	! imx(nx1:nx2,ny1:ny2,imy(nx1:nx2,ny1:ny2),ifail_L2)
40	!
41	!-----------------------------------------------------------------------
42
43	!$ USE OMP_LIB
44	use module3d_phy
45	use remplimat_declar ! les tableaux tuij, ....
46	use module_choix
47
48	! use eq_ellip_sgbsv_mod
49
50	implicit none
51
52	! declarations variables dummy
53
54	integer,intent(in) :: nx1 ! bornes du domaine en x (noeuds majeurs)
55	integer,intent(in) :: nx2
56	integer,intent(in) :: ny1 ! bornes du domaine en y (noeuds majoeurs)
57	integer,intent(in) :: ny2
58
59	!integer :: n1 ! dimension selon x
60	!integer :: n2 ! dimension selon y
61	integer :: ifail_L2 ! pour les rapports d'erreur
62
63	real, dimension(nx1:nx2,ny1:ny2),intent(inout) :: uxprec ! vitesse en entree routine
64	real, dimension(nx1:nx2,ny1:ny2),intent(inout) :: uyprec ! vitesse en entree routine
65
66	! masques vitesses.
67	! Pour donner de la souplesse dans les zones qu'on traite
68	!
69
70	! imx(i,j)=0 -> ne pas traiter cette vitesse
71	! imx(i,j)=1 -> uxnew(i,j)=uxprec(i,j) ! pour imposer les vitesses
72	! imx(i,j)=2 -> traitement général equation elliptique !
73	! imx(i,j) < 0 condition aux limites ! voir routine rempli_Tuij
74	!------------------------------------------------------------------------------------
75
76
77
78	integer, dimension(nx1:nx2,ny1:ny2),intent(inout) :: imx ! masque en entree routine
79	integer, dimension(nx1:nx2,ny1:ny2),intent(inout) :: imy ! masque en entree routine
80
81
82
83	real, dimension(nx1:nx2,ny1:ny2),intent(out) :: uxnew ! vitesse en sortie de la routine
84	real, dimension(nx1:nx2,ny1:ny2),intent(out) :: uynew ! vitesse en sortie de la routine
85
86
87	! variables locales.
88	!--------------------
89
90	real :: dx2=dx*dx ! variable de travail
91	real :: beta ! pour le frottement
92	real :: scal ! pour le conditionnement (diagonale=1)
93
94	character(len=80) :: fileout ! pour sortie fichier verification
95
96
97	! pour les fronts, on suppose que l'epaisseur est celle du noeud amont
98	! dans opposx et opposy
99
100	real, dimension(nx,ny) :: Hmx_oppos !
101	real, dimension(nx,ny) :: Hmy_oppos !
102
103	if (itracebug.eq.1) call tracebug(' Subroutine rempli_L2')
104
105
106	!-----------------------------
107	!$OMP PARALLEL
108	!$OMP WORKSHARE
109	Tu(:,:,:,:) = 0. ; Tv(:,:,:,:) = 0. ; Su(:,:,:,:) = 0. ; Sv(:,:,:,:) = 0.
110	opposx(:,:) = 0. ; opposy(:,:) = 0.
111	Mu(:,:,:,:) = 0. ; Mv(:,:,:,:) = 0. ; Nu(:,:,:,:) = 0. ; Nv(:,:,:,:) = 0.
112
113	ligu_L2(:,:) = 0 ; ligv_L2(:,:) = 0
114
115	ok_umat(:,:) = .true. ; ok_vmat(:,:) = .true.
116	ghost_x(:,:) = .false. ; ghost_y(:,:) = .false.
117
118	pos_ligu(:,:)=-9999 ; pos_ligv(:,:)=-9999
119
120	Hmx_oppos(:,:) = Hmx(:,:)
121	Hmy_oppos(:,:) = Hmy(:,:)
122	!$OMP END WORKSHARE
123	!$OMP END PARALLEL
124
125	! calcul de Hmx_oppos et Hmy_oppos dans le cas de fronts
126	! ce calcul pourrait être dans diagno et faire passer hmx_oppos par module
127
128	! 4 avril 2012 : Finalement il me semble que c'est faux, il suffit de garder le Hmx habituel.
129	! cat
130
131	!!$do j=ny1,ny2
132	!!$ do i=nx1,nx2
133	!!$
134	!!$ il=max(1,i-1)
135	!!$ jl=max(1,j-1)
136	!!$
137	!!$
138	!!$ if ((H(i,j).gt.1.).and.(H(il,j).le.1.)) then ! Bord West : on garde H(i,j)
139	!!$ Hmx_oppos(i,j)=H(i,j)
140	!!$
141	!!$ else if ((H(i,j).le.1.).and.(H(il,j).gt.1.)) then ! Bord Est : on garde H(il,j)
142	!!$ Hmx_oppos(i,j)=H(il,j)
143	!!$
144	!!$ end if
145	!!$
146	!!$ if ((H(i,j).gt.1.).and.(H(i,jl).le.1.)) then ! Bord Sud : on garde H(i,j)
147	!!$ Hmy_oppos(i,j)=H(i,j)
148	!!$
149	!!$ else if ((H(i,j).le.1.).and.(H(i,jl).gt.1.)) then ! Bord Nord on garde H(i,jl)
150	!!$ Hmy_oppos(i,j)=H(i,jl)
151	!!$
152	!!$ end if
153	!!$
154	!!$
155	!!$ end do
156	!!$ end do
157
158
159	! mailles mineures en dehors du domaine
160	!imy(:,ny1) = 0
161	!imx(nx1,:) = 0
162
163
164	! limite la valeur de frotmx a betamax
165	! devrait plutot etre fait dans dragging
166	!------------------------------------------
167	if (iter_beta.eq.0) then
168	call limit_frotm
169	else
170	!$OMP PARALLEL
171	!$OMP WORKSHARE
172	frotmx(:,:)= betamx(:,:)
173	frotmy(:,:)= betamy(:,:)
174	!$OMP END WORKSHARE
175	!$OMP END PARALLEL
176	end if
177
178
179
180	! remplissage des Tuij, ....
181	!----------------------------
182	call rempli_Tuij
183
184	! debug : impression des Tuij en un point
185	!write(166,*)
186	!write(166,*) 'impression des Tuij'
187
188	!j = 12
189
190	!do i = 244,245
191
192	!!$do i = 244,246
193	!!$ write(166,801) i,imx(i,j),frotmx(i,j),hmx_oppos(i,j),sdx(i,j)
194	!!$ do k= 1,-1,-1
195	!!$ write(166,800) k,Tu(i,j,-1:1,k),opposx(i,j)
196	!!$ end do
197	!!$ write(166,*)
198	!!$ do k= 1,-1,-1
199	!!$ write(166,800) k,Tv(i,j,-1:1,k)
200	!!$ end do
201	!!$ write(166,*)
202	!!$ write(166,*) 'impression des Svij'
203	!!$ write(166,801) i,imy(i,j),frotmy(i,j),hmy_oppos(i,j),sdy(i,j)
204	!!$ do k= 1,-1,-1
205	!!$ write(166,800) k,Su(i,j,-1:1,k),opposy(i,j)
206	!!$ end do
207	!!$ write(166,*)
208	!!$ do k= 1,-1,-1
209	!!$ write(166,800) k,Sv(i,j,-1:1,k)
210	!!$ end do
211	!!$
212	!!$
213	!!$
214	!!$800 format(' k=',(i0,1x),': ',4(es12.4,1x) )
215	!!$801 format(' i,im,frotm,hm,sd ',2(i0.1,x),3(es12.4,1x))
216	!!$ write(166,*)
217	!!$end do
218	!!$write(166,*)
219
220
221	! debug : differents termes de l equation MacAyeal
222	!---------------------------------------------------
223
224
225	!!$do j=ny1+1,ny2-1
226	!!$ do i=nx1+1,nx2-1
227	!!$
228	!!$ if (gzmx(i,j)) then
229	!!$ do jl=-1,1
230	!!$ do il=-1,1
231	!!$ debug_3D(i,j,82) =debug_3D(i,j,82)+(Tu(i,j,il,jl)*Vcol_x(i+il,j+jl) &
232	!!$ +Tv(i,j,il,jl)*Vcol_y(i+il,j+jl))
233	!!$
234	!!$ end do
235	!!$ end do
236	!!$! debug_3D(i,j,83)=frotmx(i,j)dx2Vcol_x(i,j)
237	!!$ end if
238	!!$ end do
239	!!$end do
240
241
242	if (iter_beta.eq.1) then
243	call calc_beta(Vcol_x, Vcol_y)
244	fileout='beta-estime.dat'
245	call write_datfile3(nx,ny,betamx,betamy,beta_centre,fileout)
246	if (itracebug.eq.1) call tracebug(' remplimat apres calc_beta')
247	iter_beta = 0 ! pour faire un calcul de vitesse le coup d'apres
248	corr_def = .true.
249	goto 888
250
251	else if ((iter_beta.eq.0).and.(corr_def)) then
252	call calc_beta(Vsl_x, Vsl_y)
253	endif
254
255	if (itracebug.eq.1) write(num_tracebug,*) 'iter_beta, corr_deff',iter_beta, corr_def
256
257	!!$do j=ny1,ny2
258	!!$ do i=nx1,nx2
259	!!$ debug_3D(i,j,84) = imx(i,j)
260	!!$ end do
261	!!$end do
262
263	! call rempli_Tuij ! pour reutiliser le nouveau beta
264
265
266	! Pour determiner les colonnes Mu, Mv, Nu, Nv
267	!---------------------------------------------
268	call Mu_Mv
269
270
271	! pour faire une sortie graphique de la matrice L2 et bande couchee avant tout rangement
272	!---------------------------------------------------------------------------------------
273	!call graphique_L2(2ny+2,2ny+2,nx1,nx2,ny1,ny2,imx(nx1:nx2,ny1:ny2),imy(nx1:nx2,ny1:ny2))
274
275
276	! pour determiner les noeuds qui participent a l'equation elliptique
277	!-------------------------------------------------------------------
278	call okmat0
279
280	! pour avoir la diagonale=1 a faire imperativement APRES okmat
281	!--------------------------------------------------------------
282	!$OMP PARALLEL
283	!$OMP DO PRIVATE(scal)
284	do j=ny1,ny2
285	do i=nx1,nx2
286
287
288	if (imx(i,j).ne.0) then ! test pour eviter les noeuds non traites
289
290	scal=Tu(i,j,0,0)
291	if(abs(scal).lt.1.e-30) then
292	write(num_tracebug,*)' i,j,imx,pvi',i,j,imx(i,j),pvi(i,j)
293	! write(num_tracebug,*) Tu(i,j,:,:)
294	end if
295
296
297	Tu(i,j,:,:)=Tu(i,j,:,:)/scal ! qui produiraient une division par 0
298	Tv(i,j,:,:)=Tv(i,j,:,:)/scal ! Tu(:,:,0,0) est la diagonale
299	opposx(i,j)=opposx(i,j)/scal
300
301	! if ((i.eq.200).and.(j.eq.251)) then
302	! write(num_tracebug,*)' i,j,imx',i,j,imx(i,j)
303	! write(num_tracebug,*)' pvi,scal',pvi(i,j),scal
304	! write(num_tracebug,*) Tu(i,j,:,-1)
305	! write(num_tracebug,*) Tu(i,j,:,0)
306	! write(num_tracebug,*) Tu(i,j,:,1)
307	! write(num_tracebug,*)' opposx',opposx(i,j)
308	! end if
309
310
311
312	end if
313
314	if (imy(i,j).ne.0) then ! test pour eviter les noeuds non traites
315	scal=Sv(i,j,0,0)
316	if(abs(scal).lt.1.e-30) then
317	write(num_tracebug,*)' i,j,Svij,imy,pvi',i,j,imy(i,j),pvi(i,j)
318	! write(num_tracebug,*) Sv(i,j,:,:)
319	end if
320
321	Su(i,j,:,:)=Su(i,j,:,:)/scal ! qui produiraient une division par 0
322	Sv(i,j,:,:)=Sv(i,j,:,:)/scal ! Sv(:,:,0,0) est la diagonale
323	opposy(i,j)=opposy(i,j)/scal
324	end if
325
326
327
328	end do
329	end do
330	!$OMP END DO
331	!$OMP END PARALLEL
332
333	! mise a identite des noeuds fantomes
334	call ghost_identite
335
336	call ghost_remove
337
338	! call graphique_L2(kl,ku,nx1,nx2,ny1,ny2,imx(nx1:nx2,ny1:ny2),imy(nx1:nx2,ny1:ny2))
339
340
341	!debug_3d(:,:,35) = opposx(:,:)
342	!debug_3d(:,:,36) = opposy(:,:)
343
344
345
346	call resol_ellipt(nx1,nx2,ny1,ny2, & ! bornes du domaine
347
348	uxprec(nx1:nx2,ny1:ny2), uyprec(nx1:nx2,ny1:ny2), & ! vitesses precedentes
349
350
351	uxnew(nx1:nx2,ny1:ny2), uynew(nx1:nx2,ny1:ny2), & ! nouvelles vitesses
352
353	imx(nx1:nx2,ny1:ny2), imy(nx1:nx2,ny1:ny2),ifail_L2) ! masques
354
355	if (itracebug.eq.1) call tracebug ('apres subroutine resol_ellipt')
356
357	if (itracebug.eq.1) write(num_tracebug,*) '(ifail_L2 (si 0 pas de probleme)', ifail_L2
358
359	! remise a 0 des noeuds fantomes
360	!$OMP PARALLEL
361	!$OMP WORKSHARE
362	where ((Hmx(nx1:nx2,ny1:ny2).le.1.).and.(flgzmx(nx1:nx2,ny1:ny2)))
363	uxnew(nx1:nx2,ny1:ny2)=0.
364	end where
365
366	where ((Hmy(nx1:nx2,ny1:ny2).le.1.).and.(flgzmy(nx1:nx2,ny1:ny2)))
367	uynew(nx1:nx2,ny1:ny2)=0.
368	end where
369	!$OMP END WORKSHARE
370	!$OMP END PARALLEL
371
372	! call graphique_L2(kl,ku,nx1,nx2,ny1,ny2,imx(nx1:nx2,ny1:ny2),imy(nx1:nx2,ny1:ny2))
373
374
375
376	! pour appeler la routine de resolution du systeme lineaire
377
378	!!$i = 244
379	!!$do j=11,11,-1
380	!!$
381	!!$write(166,*) 'j=',j
382	!!$ write(166,803) i,uxprec(i,j),i+1,uxprec(i+1,j),i+2,uxprec(i+2,j)
383	!!$ write(166,802) i,uxnew(i,j),i+1,uxnew(i+1,j),i+2,uxnew(i+2,j)
384	!!$802 format(3('uxnew ',' i=',(i0,1x,es12.4,1x) ))
385	!!$803 format(3('uxprec',' i=',(i0,1x,es12.4,1x) ))
386	!!$write(166,*)
387
388	!!$end do
389
390	! limitation des nouvelles vitesses
391
392	!!$uxnew(:,:)=max(-3900.,uxnew(:,:))
393	!!$uxnew(:,:)=min( 3900.,uxnew(:,:))
394	!!$uynew(:,:)=max(-3900.,uynew(:,:))
395	!!$uynew(:,:)=min( 3900.,uynew(:,:))
396	!!$
397
398
399
400
401	888 return
402
403
404
405
406	contains ! pour que les subroutines qui suivent partagent les variables de rempli_L2
407
408
409	!-----------------------------------------------------------------------
410	! routines internes au module
411	!
412	! tout ce qui suite est ecrit sous forme de routines
413	! pour rendre plus lisible remplidom
414
415
416	subroutine limit_frotm
417
418	!$USE OMP_LIB
419	!-------------------------
420	! limite la valeur de frotmx a betamax
421	! devrait plutot etre fait dans dragging
422
423	if (itracebug.eq.1) call tracebug(' Subroutine limit_frotm')
424	if (itracebug.eq.1) write(Num_tracebug,*)'betamax = ',betamax
425
426	!$OMP PARALLEL
427	!$OMP WORKSHARE
428	where (flgzmx(:,:))
429	frotmx(:,:)=min(abs(betamx(:,:)),betamax_2d(:,:))
430	elsewhere
431	frotmx(:,:)=0
432	end where
433
434	where (flgzmy(:,:))
435	frotmy(:,:)=min(abs(betamy(:,:)),betamax_2d(:,:))
436	elsewhere
437	frotmy(:,:)=0
438	end where
439	!$OMP END WORKSHARE
440	!$OMP END PARALLEL
441
442	end subroutine limit_frotm
443	!------------------------------------------------------------------
444
445
446	!------------------------------------------------------------------
447	!
448	! subroutines calcul des Tuij Cat juin 2008
449	!
450
451	subroutine rempli_Tuij ! appelle tous les cas
452
453	!$ USE OMP_LIB
454
455	! imx(i,j)=0 -> ne pas traiter cette vitesse
456	! imx(i,j)=1 -> uxnew(i,j)=uxprec(i,j) ! pour imposer les vitesses
457	! imx(i,j)=2 -> traitement général equation elliptique !
458
459	! imx(i,j) < 0 condition aux limites
460	!------------------------------------
461	! numerotation en tournant dans le sens direct.
462	! -1 bord Sud (en bas), -2 bord Est (a droite)
463	! -3 bord Nord (en haut), -4 bord West (a gauche).
464	! les coins ont l'indice combine -12 coin en bas a gauche
465
466	! -34 -3 Nord -23
467	! !----------------------------------------!
468	! ! !
469	! ! 1 (prescrite) !
470	! -4 ! ou ! -2
471	! West! 2 (L2) ! Est
472	! ! !
473	! ! !
474	! !----------------------------------------!
475	! -41 -1 Sud -12
476
477
478	if (itracebug.eq.1) call tracebug(' Rempli Tuij')
479
480
481
482
483	count_line=1 ! pour numeroter les lignes
484	!------------------------------------------------------------------
485	!$OMP PARALLEL
486	!$OMP DO
487	lignes_UV: do j=ny1,ny2
488	do i=nx1,nx2
489
490	if (i.gt.nx1) then ! Vitesses U
491	ligu_L2(i,j)=count_line
492	!$OMP ATOMIC
493	count_line=count_line+1
494
495	case_imx: select case (imx(i,j)) ! la routine appelee depend d'imx
496	!------------------------------------
497
498	case(0) ! ne pas traiter ce point
499	ligu_L2(i,j)=0
500	count_line=count_line-1
501
502	case(1) ! vitesse imposee
503	call vel_U_presc(uxprec(i,j),Tu(i,j,0,0),opposx(i,j))
504
505	case(2) ! cas general
506	call vel_U_general(frotmx(i,j),dx2,pvi(i,j),pvi(i-1,j),pvm(i,j),pvm(i,j+1),sdx(i,j),rog,hmx_oppos(i,j), &
507	Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
508
509	case(-1) ! bord sud
510	call vel_U_South(Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
511
512	case(-2) ! bord Est
513	call vel_U_East(pvi(i,j),pvi(i-1,j),rog,H(i-1,j),rowg,slv(i-1,j),B(i-1,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
514
515	case(-3) ! bord nord
516	call vel_U_North(Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
517
518	case(-4) ! bord West
519	call vel_U_West(pvi(i,j),pvi(i-1,j),rog,H(i,j),rowg,slv(i,j),B(i,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
520
521	case(-12) ! coin SE
522	call vel_U_SE(Tu(i,j,:,:),opposx(i,j))
523
524	case(-23) ! coin NE
525	call vel_U_NE(Tu(i,j,:,:),opposx(i,j))
526
527	case(-34) ! coin NW
528	call vel_U_NW(Tu(i,j,:,:),opposx(i,j))
529
530	case(-41) ! coin SW
531	call vel_U_SW(Tu(i,j,:,:),opposx(i,j))
532
533	end select case_imx
534	end if
535
536	if (j.gt.ny1) then ! Vitesses V
537
538	ligv_L2(i,j)=count_line
539	!$OMP ATOMIC
540	count_line=count_line+1
541
542	case_imy: select case (imy(i,j)) ! la routine appelee depend d'imy
543	!------------------------------------
544
545	case(0) ! ne pas traiter ce point
546	ligv_L2(i,j)=0
547	count_line=count_line-1
548
549	case(1) ! vitesse imposee
550	call vel_V_presc(uyprec(i,j),Sv(i,j,0,0),opposy(i,j))
551
552	case(2) ! cas general
553	call vel_V_general(frotmy(i,j),dx2,pvi(i,j),pvi(i,j-1),pvm(i,j),pvm(i+1,j),sdy(i,j),rog,hmy_oppos(i,j), &
554	Su(i,j,:,:),Sv(i,j,:,:),opposy(i,j))
555
556	case(-1) ! bord sud
557	call vel_V_South(pvi(i,j),pvi(i,j-1),rog,H(i,j),rowg,slv(i,j),B(i,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
558
559	case(-2) ! bord Est
560	call vel_V_East(Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
561
562	case(-3) ! bord nord
563	call vel_V_North(pvi(i,j),pvi(i,j-1),rog,H(i-1,j),rowg,slv(i-1,j),B(i-1,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
564
565	case(-4) ! bord West
566	call vel_V_West(Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
567
568	case(-12) ! coin SE
569	call vel_V_SE(Sv(i,j,:,:),opposy(i,j))
570
571	case(-23) ! coin NE
572	call vel_V_NE(Sv(i,j,:,:),opposy(i,j))
573
574	case(-34) ! coin NW
575	call vel_V_NW(Sv(i,j,:,:),opposy(i,j))
576
577	case(-41) ! coin SW
578
579	call vel_V_SW(Sv(i,j,:,:),opposy(i,j))
580
581	end select case_imy
582	end if
583	end do
584	end do lignes_UV
585	!$OMP END DO
586	!$OMP END PARALLEL
587
588	end subroutine rempli_Tuij
589	!------------------------------------------------------------------
590
591	! Les differents cas :
592	! Vitesses U, puis vitesses V les coins sont traites a la fin
593
594
595	! vitesses U
596	!------------------------------------------------------------------
597
598	subroutine vel_U_presc(uxprec,Tu,opposx) ! ! vitesse imposee
599
600	implicit none
601
602	real,intent(in) :: uxprec
603	real,intent(out) :: Tu
604	real,intent(out) :: opposx
605
606	Tu=1.
607	opposx=uxprec
608
609	end subroutine vel_U_presc
610	!------------------------------------------------------------------
611
612	subroutine vel_U_general(frotmx,dx2,pvi,pvi_imoinsun,pvm,pvm_jplusun,sdx,rog,hmx_oppos,Tu,Tv,opposx) ! cas general
613
614	implicit none
615
616	real,intent(in) :: frotmx
617	real,intent(in) :: dx2
618	real,intent(in) :: pvi
619	real,intent(in) :: pvi_imoinsun
620	real,intent(in) :: pvm
621	real,intent(in) :: pvm_jplusun
622	real,intent(in) :: sdx
623	real,intent(in) :: rog
624	real,intent(in) :: hmx_oppos
625	real,dimension(-1:1,-1:1),intent(inout) :: Tu
626	real,dimension(-1:1,-1:1),intent(inout) :: Tv
627	real,intent(out) :: opposx
628
629	beta=frotmx*dx2 ! Terme en u(i,j)
630
631	Tu(0,0) = -4.pvi - 4.pvi_imoinsun - pvm_jplusun - pvm - beta
632
633	Tu(-1,0) = 4.*pvi_imoinsun ! Terme en u(i-1,j)
634
635	Tu(1,0) = 4.*pvi ! Terme en u(i+1,j)
636
637	Tu(0,-1) = pvm ! Terme en u(i,j-1)
638
639	Tu(0,1) = pvm_jplusun ! Terme en u(i,j+1)
640
641	Tv(0,0) = -2.*pvi - pvm ! Terme en v(i,j)
642
643	Tv(-1,0) = 2.*pvi_imoinsun + pvm ! Terme en v(i-1,j)
644
645	Tv(-1,1) = -2.*pvi_imoinsun - pvm_jplusun ! Terme en v(i-1,j+1)
646
647	Tv(0,1) = 2.*pvi + pvm_jplusun ! Terme en v(i,j+1)
648
649	opposx = roghmx_oppossdx*dx2 ! Terme en opposx(i,j)
650
651	return
652	end subroutine vel_U_general
653
654	!------------------------------------------------------------------
655	! voir explications dans vel_U_West
656
657	subroutine vel_U_South(Tu,Tv,opposx) ! bord sud non cisaillement
658
659	implicit none
660	real,dimension(-1:1,-1:1),intent(inout) :: Tu
661	real,dimension(-1:1,-1:1),intent(inout) :: Tv
662	real, intent(out) :: opposx
663
664	Tu(0,0) = 1.
665	Tu(0,1) = -1.
666	Tv(0,1) = -1.
667	Tv(-1,1) = 1.
668	opposx = 0.
669
670	return
671	end subroutine vel_U_South
672	!------------------------------------------------------------------
673
674	subroutine vel_U_North(Tu,Tv,opposx) ! bord nord non cisaillement
675
676	implicit none
677	real,dimension(-1:1,-1:1),intent(inout) :: Tu
678	real,dimension(-1:1,-1:1),intent(inout) :: Tv
679	real, intent(out) :: opposx
680
681	Tu(0,0) = 1.
682	Tu(0,-1) = -1.
683	Tv(0,0) = 1.
684	Tv(-1,0) = -1.
685	opposx = 0.
686
687	return
688	end subroutine vel_U_North
689	!------------------------------------------------------------------
690	! voir explications dans vel_U_West
691
692	subroutine vel_U_East(pvi,pvi_imoinsun,rog,H_imoinsun,rowg,slv_imoinsun,B_imoinsun,dx,Tu,Tv,opposx) ! bord Est pression
693
694	implicit none
695	real,intent(in) :: pvi
696	real,intent(in) :: pvi_imoinsun
697	real,intent(in) :: rog
698	real,intent(in) :: H_imoinsun
699	real,intent(in) :: rowg
700	real,intent(in) :: slv_imoinsun
701	real,intent(in) :: B_imoinsun
702	real,intent(in) :: dx
703	real,dimension(-1:1,-1:1),intent(inout) :: Tu
704	real,dimension(-1:1,-1:1),intent(inout) :: Tv
705	real, intent(out) :: opposx
706
707	Tu(0,0) = 4.*pvi_imoinsun
708	Tu(-1,0) = -4.*pvi_imoinsun
709	Tv(0,1) = 2.*pvi
710	Tv(0,0) = -2.*pvi
711	opposx = 0.5(rogH_imoinsun*2-rowg(max(slv_imoinsun-B_imoinsun,0.))*2)dx
712
713	return
714	end subroutine vel_U_East
715	!------------------------------------------------------------------
716
717	subroutine vel_U_West(pvi,pvi_imoinsun,rog,H,rowg,slv,B,dx,Tu,Tv,opposx) ! bord West pression
718	! tous les coef * -1
719	implicit none
720	real,intent(in) :: pvi
721	real,intent(in) :: pvi_imoinsun
722	real,intent(in) :: rog
723	real,intent(in) :: H
724	real,intent(in) :: rowg
725	real,intent(in) :: slv
726	real,intent(in) :: B
727	real,intent(in) :: dx
728	real,dimension(-1:1,-1:1),intent(inout) :: Tu
729	real,dimension(-1:1,-1:1),intent(inout) :: Tv
730	real, intent(out) :: opposx
731
732
733	! Tu(i,j,0,0) = 4.*pvi(i-1,j)
734	! Tu(i,j,1,0) = -4.*pvi(i-1,j)
735
736	Tu(0,0) = 4.*pvi ! quand le bord est epais, il faut prendre
737	Tu(1,0) = -4.*pvi ! le noeud du shelf pas celui en H=1
738	! s'il est fin, pas de difference
739
740	!!$ Tv(i-1,j,-1,0) = 2.*pvi(i-1,j) ! ????
741	!!$ Tv(i-1,j,-1,1) = -2.*pvi(i-1,j)
742
743	Tv(-1,0) = 2.*pvi_imoinsun
744	Tv(-1,1) = -2.*pvi_imoinsun
745	opposx = -0.5(rogH*2-rowg(max(slv-B,0.))*2)dx
746
747
748	return
749	end subroutine vel_U_West
750	!------------------------------------------------------------------
751
752
753	! vitesses V
754	!------------------------------------------------------------------
755
756	subroutine vel_V_presc(uyprec,Sv,opposy) ! vitesse imposee
757
758	implicit none
759	real,intent(in) :: uyprec
760	real,intent(out) :: Sv
761	real,intent(out) :: opposy
762
763	Sv=1.
764	opposy=uyprec
765
766	return
767	end subroutine vel_V_presc
768	!------------------------------------------------------------------
769
770	subroutine vel_V_general(frotmy,dx2,pvi,pvi_jmoinsun,pvm,pvm_iplusun,sdy,rog,hmy_oppos,Su,Sv,opposy)! cas general
771
772	implicit none
773	real,intent(in) :: frotmy
774	real,intent(in) :: dx2
775	real,intent(in) :: pvi
776	real,intent(in) :: pvi_jmoinsun
777	real,intent(in) :: pvm
778	real,intent(in) :: pvm_iplusun
779	real,intent(in) :: sdy
780	real,intent(in) :: rog
781	real,intent(in) :: hmy_oppos
782	real,dimension(-1:1,-1:1),intent(inout) :: Su
783	real,dimension(-1:1,-1:1),intent(inout) :: Sv
784	real,intent(out) :: opposy
785
786	real :: beta
787
788
789	beta = frotmy*dx2 ! Terme en v(i,j)
790	Sv(0,0) = -4.pvi - 4.pvi_jmoinsun - pvm_iplusun - pvm-beta
791
792	Sv(0,-1) = 4.*pvi_jmoinsun ! Terme en v(i,j-1)
793
794	Sv(0,1) = 4.*pvi ! Terme en v(i,j+1)
795
796	Sv(-1,0) = pvm ! Terme en v(i-1,j)
797
798	Sv(1,0) = pvm_iplusun ! Terme en v(i+1,j)
799
800	Su(0,0) = -2.*pvi - pvm ! Terme en u(i,j)
801
802	Su(0,-1) = 2.*pvi_jmoinsun + pvm ! Terme en u(i,j-1)
803
804	Su(1,-1) = -2.*pvi_jmoinsun - pvm_iplusun ! Terme en u(i+1,j-1)
805
806	Su(1,0) = 2.*pvi + pvm_iplusun ! Terme en u(i+1,j)
807
808	opposy = roghmy_oppossdy*dx2 ! Terme en opposy(i,j)
809
810	return
811	end subroutine vel_V_general
812
813	!------------------------------------------------------------------
814
815	subroutine vel_V_West(Sv,Su,opposy) ! bord West non cisaillement
816
817	implicit none
818	real,dimension(-1:1,-1:1),intent(inout) :: Sv
819	real,dimension(-1:1,-1:1),intent(inout) :: Su
820	real,intent(out) :: opposy
821
822	Sv(0,0) = 1.
823	Sv(1,0) = -1.
824	Su(1,0) = -1.
825	Su(1,-1) = 1.
826	opposy = 0.
827
828	return
829	end subroutine vel_V_West
830	!------------------------------------------------------------------
831
832
833	subroutine vel_V_East(Sv,Su,opposy) ! bord Est non cisaillement
834
835	implicit none
836	real,dimension(-1:1,-1:1),intent(inout) :: Sv
837	real,dimension(-1:1,-1:1),intent(inout) :: Su
838	real,intent(out) :: opposy
839
840	Sv(0,0) = 1.
841	Sv(-1,0) = -1.
842	Su(0,0) = 1.
843	Su(0,-1) = -1.
844	opposy = 0.
845
846	return
847	end subroutine vel_V_East
848	!------------------------------------------------------------------
849	!-----------------------------------------------------------------------
850	! voir explications dans vel_U_West
851
852	subroutine vel_V_North(pvi,pvi_jmoinsun,rog,H_imoinsun,rowg,slv_imoinsun,B_imoinsun,dx,Sv,Su,opposy) ! bord Nord pression
853
854	implicit none
855	real,intent(in) :: pvi
856	real,intent(in) :: pvi_jmoinsun
857	real,intent(in) :: rog
858	real,intent(in) :: H_imoinsun
859	real,intent(in) :: rowg
860	real,intent(in) :: slv_imoinsun
861	real,intent(in) :: B_imoinsun
862	real,intent(in) :: dx
863	real,dimension(-1:1,-1:1),intent(inout) :: Sv
864	real,dimension(-1:1,-1:1),intent(inout) :: Su
865	real,intent(out) :: opposy
866
867	Sv(0,0) = 4.*pvi_jmoinsun
868	Sv(0,-1) = -4.*pvi_jmoinsun
869	Su(1,0) = 2.*pvi
870	Su(0,0) = -2.*pvi
871	opposy = 0.5(rogH_imoinsun*2-rowg(max(slv_imoinsun-B_imoinsun,0.))*2)dx
872
873	return
874	end subroutine vel_V_North
875	!------------------------------------------------------------------
876	! voir explications dans vel_U_West
877
878	subroutine vel_V_South(pvi,pvi_jmoinsun,rog,H,rowg,slv,B,dx,Sv,Su,opposy) ! bord sud pression
879	! tous les coef * -1
880	implicit none
881	real,intent(in) :: pvi
882	real,intent(in) :: pvi_jmoinsun
883	real,intent(in) :: rog
884	real,intent(in) :: H
885	real,intent(in) :: rowg
886	real,intent(in) :: slv
887	real,intent(in) :: B
888	real,intent(in) :: dx
889	real,dimension(-1:1,-1:1),intent(inout) :: Sv
890	real,dimension(-1:1,-1:1),intent(inout) :: Su
891	real,intent(out) :: opposy
892
893	Sv(0,0) = 4.*pvi
894	Sv(0,1) =-4.*pvi
895	Su(0,-1) = 2.*pvi_jmoinsun
896	Su(1,-1) =-2.*pvi_jmoinsun
897	opposy = -0.5(rogH*2-rowg(max(slv-B,0.))*2)dx
898
899	return
900	end subroutine vel_V_South
901	!------------------------------------------------------------------
902
903
904	! coins
905	!------------------------------------------------------------------
906
907	! Dans les coins, la condition non-cisaillement est la meme
908	! pour les noeuds en U et en V. Il ne faut donc l'utiliser qu'une seule fois
909	! On ajoute alors comme condition la symetrie
910	!
911	! du coup la condition s'ecrit dU/dy=0 en U et dV/dx=0 en V
912	! la condition de non cisaillement est automatiquement satisfaite
913	!
914
915
916	!Coin SE
917	!----------
918	subroutine vel_U_SE(Tu,opposx)
919
920	implicit none
921	real,dimension(-1:1,-1:1),intent(inout) :: Tu
922	real,intent(out) :: opposx
923
924	Tu(0,0) = 1. !*pvi(i,j)
925	Tu(0,1) = -1. !*pvi(i,j)
926	opposx = 0.
927
928	return
929	end subroutine vel_U_SE
930	!------------------------------------------------------------------
931
932	subroutine vel_V_SE(Sv,opposy)
933
934	implicit none
935	real,dimension(-1:1,-1:1),intent(inout) :: Sv
936	real,intent(out) :: opposy
937
938	Sv(0,0) = 1. !*pvi(i,j)
939	Sv(-1,0)= -1. !*pvi(i,j)
940	opposy = 0.
941
942	return
943	end subroutine vel_V_SE
944	!------------------------------------------------------------------
945
946
947	!Coin SW
948	!----------
949	subroutine vel_U_SW(Tu,opposx)
950
951	implicit none
952	real,dimension(-1:1,-1:1),intent(inout) :: Tu
953	real,intent(out) :: opposx
954
955	Tu(0,0) = 1. ! *pvi(i,j)
956	Tu(0,1) = -1. ! *pvi(i,j)
957	opposx = 0.
958
959	return
960	end subroutine vel_U_SW
961	!------------------------------------------------------------------
962
963	subroutine vel_V_SW(Sv,opposy)
964
965	implicit none
966	real,dimension(-1:1,-1:1),intent(inout) :: Sv
967	real,intent(out) :: opposy
968
969	Sv(0,0) = 1. !*pvi(i,j)
970	Sv(1,0) = -1. !*pvi(i,j)
971	opposy = 0.
972
973	return
974	end subroutine vel_V_SW
975	!------------------------------------------------------------------
976
977	!Coin NE
978	!----------
979	subroutine vel_U_NE(Tu,opposx)
980
981	implicit none
982	real,dimension(-1:1,-1:1),intent(inout) :: Tu
983	real,intent(out) :: opposx
984
985	Tu(0,0) = 1. ! *pvi(i,j)
986	Tu(0,-1) = -1. ! *pvi(i,j)
987	opposx = 0.
988
989	return
990	end subroutine vel_U_NE
991	!------------------------------------------------------------------
992
993	subroutine vel_V_NE(Sv,opposy)
994
995	implicit none
996	real,dimension(-1:1,-1:1),intent(inout) :: Sv
997	real,intent(out) :: opposy
998
999	Sv(0,0) = 1. ! *pvi(i,j)
1000	Sv(-1,0) = -1. ! *pvi(i,j)
1001	opposy = 0.
1002
1003	return
1004	end subroutine vel_V_NE
1005	!------------------------------------------------------------------
1006
1007	!Coin NW
1008	!----------
1009	subroutine vel_U_NW(Tu,opposx)
1010
1011	implicit none
1012	real,dimension(-1:1,-1:1),intent(inout) :: Tu
1013	real,intent(out) :: opposx
1014
1015	Tu(0,0) = 1. ! *pvi(i,j)
1016	Tu(0,-1) = -1. ! *pvi(i,j)
1017	opposx = 0.
1018
1019	return
1020	end subroutine vel_U_NW
1021	!------------------------------------------------------------------
1022
1023	subroutine vel_V_NW(Sv,opposy)
1024
1025	implicit none
1026	real,dimension(-1:1,-1:1),intent(inout) :: Sv
1027	real,intent(out) :: opposy
1028
1029	Sv(0,0) = 1. !*pvi(i,j)
1030	Sv(1,0) = -1. !*pvi(i,j)
1031	opposy = 0.
1032
1033	return
1034	end subroutine vel_V_NW
1035	!------------------------------------------------------------------
1036
1037
1038	! fin des routines Tuij,...
1039	!------------------------------------------------------------------
1040
1041
1042
1043	subroutine Mu_Mv
1044
1045	!---------------------------------------------------------
1046	! Pour déterminer Mu, Mv, Nu, Nv les colonnes de L2
1047	! A appeler avant la division des Tu, Sv par la diagonale
1048	!---------------------------------------------------------
1049
1050	! pourrait être dans un autre fichier
1051
1052	if (itracebug.eq.1) call tracebug(' Mu_Mv')
1053
1054	Col_U: do j=ny1,ny2 ! balaye tous les noeuds U
1055	do i=nx1+1,nx2
1056
1057	! ilmin et jlmin boucle U
1058	ilmin=max(-1,nx1+1-i) ! pour avoir i+il entre nx1+1 et nx2
1059	ilmax=min(1,nx2-i)
1060
1061	jlmin=max(-1,ny1-j) ! pour avoir j+jl entre ny1 et ny2
1062	jlmax=min(1,ny2-j)
1063
1064
1065
1066	do jl=jlmin,jlmax ! balaye tous les coefficients d'un noeud U
1067	do il=ilmin,ilmax
1068
1069	Mu(i+il,j+jl,-il,-jl)=abs(Tu(i,j,il,jl)) ! abs pour faciliter les tests.
1070	Mv(i+il,j+jl,-il,-jl)=abs(Tv(i,j,il,jl))
1071
1072
1073	end do
1074	end do
1075	diagU(i,j)=Tu(i,j,0,0)
1076
1077	end do
1078	end do Col_U
1079
1080	Col_V: do j=ny1+1,ny2 ! balaye tous les noeuds V
1081	do i=nx1,nx2
1082
1083	! ilmin et jlmin boucle V
1084	ilmin=max(-1,nx1-i) ! pour avoir i+il entre nx1 et nx2
1085	ilmax=min(1,nx2-i)
1086
1087	jlmin=max(-1,ny1+1-j) ! pour avoir j+jl entre ny1+1 et ny2
1088	jlmax=min(1,ny2-j)
1089
1090	do jl=jlmin,jlmax ! balaye tous les coefficients d'un noeud U
1091	do il=ilmin,ilmax
1092
1093	Nu(i+il,j+jl,-il,-jl)=abs(Su(i,j,il,jl)) ! abs pour faciliter les tests.
1094	Nv(i+il,j+jl,-il,-jl)=abs(Sv(i,j,il,jl))
1095
1096	end do
1097	end do
1098
1099	diagV(i,j)=Sv(i,j,0,0)
1100
1101
1102	end do
1103	end do Col_V
1104
1105	return
1106	end subroutine Mu_Mv
1107	!----------------------------------------------------------------------------------
1108
1109	!-------------------------------------------------------------------------
1110	subroutine okmat0
1111
1112	! pourrais être dans un autre fichier avec un appel (nx1,nx2,ny1,ny2)
1113
1114	! determine les noeuds qui participent a l'equation elliptique
1115
1116	! ce sont les noeuds dont aucun autre noeud ne dépend (colonne Mu,Nu nulle)
1117	! et qui ne dépendent d'aucun noeud (ligne Tu,Tv nulle sauf diagonale)
1118
1119
1120
1121	! tableaux de travail
1122	real, dimension(-1:1,-1:1) :: Mloc ! sous tableau local de M
1123	real, dimension(-1:1,-1:1) :: Nloc ! sous tableau local de N
1124	real, dimension(-1:1,-1:1) :: Tuloc ! sous tableau local de Tu
1125	real, dimension(-1:1,-1:1) :: Tvloc ! sous tableau local de Tv
1126	real, dimension(-1:1,-1:1) :: Suloc ! sous tableau local de Su
1127	real, dimension(-1:1,-1:1) :: Svloc ! sous tableau local de Sv
1128
1129	! test pour les noeuds fantomes
1130	real :: pvi_ghost
1131
1132
1133	if (itracebug.eq.1) call tracebug(' Okmat')
1134
1135	! initialisation
1136
1137	pvi_ghost=40.*pvimin
1138	!pvi_ghost=10.*pvimin ! condition trop faible : pas de ghost
1139
1140	ok_umat(:,:) =.false.
1141	ok_vmat(:,:) =.false.
1142	ghost_x(:,:) =.false.
1143	ghost_y(:,:) =.false.
1144
1145
1146	do j=ny1,ny2
1147	do i=nx1,nx2
1148
1149	! okumat vitesse U
1150	!---------------------------
1151
1152
1153	! condition colonne.
1154	!---------------------------
1155	Mloc(:,:)=abs(Mu(i,j,:,:))
1156	Nloc(:,:)=abs(Nu(i,j,:,:))
1157	Mloc(0,0)=0. ! diagonale
1158
1159	! okumat vrai si au moins un terme de la colonne est > eps_col et imx non 0
1160	ok_umat(i,j)=((any(Mloc.gt.eps_col)).or.(any(Nloc.gt.eps_col))).and.(imx(i,j).ne.0)
1161
1162
1163	! ghost_x est vrai si tous les elements sont inferieurs a pvi_ghost
1164	ghost_x(i,j)=((all(Mloc.lt.pvi_ghost)).and.(all(Nloc.lt.pvi_ghost)))
1165
1166
1167
1168	! condition ligne
1169	!---------------------------
1170	Tuloc(:,:)=abs(Tu(i,j,:,:))
1171	Tvloc(:,:)=abs(Tv(i,j,:,:))
1172
1173	! ghost_x est vrai si tous les elements y compris diagonale sont inferieurs a pvi_ghost
1174	ghost_x(i,j)=(ghost_x(i,j).and.((all(Tuloc.lt.pvi_ghost)).and.(all(Tvloc.lt.pvi_ghost))))
1175
1176
1177	Tuloc(0,0)=0. ! mettre cette ligne en commentaire pour avoir
1178	! aussi les noeuds identite
1179
1180	! ok_umat est vrai si au moins un des termes de la ligne est > eps_col
1181	ok_umat(i,j)=ok_umat(i,j).or.((any(Tuloc.gt.eps_col)).or.(any(Tvloc.gt.eps_col)))
1182
1183
1184
1185	! on elimine les noeuds en imx=0
1186	ok_umat(i,j)=ok_umat(i,j).and.(imx(i,j).ne.0)
1187
1188
1189	! on elimine les noeuds en ghost
1190	ok_umat(i,j)= ok_umat(i,j).and.(.not.ghost_x(i,j))
1191
1192	! okvmat vitesse V
1193	!---------------------------
1194
1195	! condition colonne.
1196	!-------------------------------
1197	Mloc(:,:)=abs(Mv(i,j,:,:))
1198	Nloc(:,:)=abs(Nv(i,j,:,:))
1199	Nloc(0,0)=0.
1200
1201	! okvmat vrai si au moins un terme de la colonne est > eps_col et imy non 0
1202	ok_vmat(i,j)=((any(Mloc.gt.eps_col)).or.(any(Nloc.gt.eps_col))).and.(imy(i,j).ne.0)
1203
1204
1205	! ghost_y est vrai si tous les elements sont inferieurs a pvi_ghost
1206	ghost_y(i,j)=((all(Mloc.lt.pvi_ghost)).and.(all(Nloc.lt.pvi_ghost)))
1207
1208
1209	! condition ligne
1210	!--------------------------------
1211	Suloc(:,:)=abs(Su(i,j,:,:))
1212	Svloc(:,:)=abs(Sv(i,j,:,:))
1213
1214	! ghost_y est vrai si tous les elements y compris diagonale sont inferieurs a pvi_ghost
1215	ghost_y(i,j)=(ghost_y(i,j).and.((all(Suloc.lt.pvi_ghost)).and.(all(Svloc.lt.pvi_ghost))))
1216
1217	Svloc(0,0)=0. ! mettre cette ligne en commentaire pour avoir
1218	! aussi les noeuds identite
1219
1220	! okvmat est vrai si au moins un des termes de la ligne est > eps_col
1221	ok_vmat(i,j)=ok_vmat(i,j).or.((any(Suloc.gt.eps_col)).or.(any(Svloc.gt.eps_col)))
1222
1223
1224	! on elimine les noeuds en imy=0
1225	ok_vmat(i,j)=ok_vmat(i,j).and.(imy(i,j).ne.0)
1226
1227	! on elimine les noeuds en ghost
1228	ok_vmat(i,j)=ok_vmat(i,j).and.(.not.ghost_y(i,j))
1229
1230	! test
1231	if ((.not.ghost_y(i,j)).and.(.not.ok_vmat(i,j))) then
1232	write(6,*)'pb ok_vmat',i,j
1233	write(6,*) 'Mloc',Mloc
1234	write(6,*) 'Nloc',Nloc
1235	write(6,*) 'Suloc',Suloc
1236	write(6,*) 'Svloc',Svloc
1237	end if
1238
1239
1240	end do
1241	end do
1242
1243
1244	! on enleve les lignes 1 decalees
1245	ok_umat(1,:)=.false.
1246	ok_vmat(:,1)=.false.
1247
1248	! sortie Netcdf pour verifier ok_umat
1249	where (ok_umat(:,:))
1250	debug_3D(:,:,39)=1
1251	elsewhere
1252	debug_3D(:,:,39)=0
1253	end where
1254
1255	where (ok_vmat(:,:))
1256	debug_3D(:,:,40)=1
1257	elsewhere
1258	debug_3D(:,:,40)=0
1259	end where
1260
1261	! sortie Netcdf pour verifier ghost
1262	where (ghost_x(:,:))
1263	debug_3D(:,:,41)=1
1264	elsewhere
1265	debug_3D(:,:,41)=0
1266
1267	end where
1268
1269	where (ghost_y(:,:))
1270	debug_3D(:,:,42)=1
1271	elsewhere
1272	debug_3D(:,:,42)=0
1273	end where
1274
1275	if (itracebug.eq.1) call tracebug(' Sortie Okmat')
1276	return
1277	end subroutine okmat0
1278
1279	!--------------------------------------------------------------------------------------
1280	subroutine ghost_identite ! mise a identite des noeuds fantomes
1281
1282	do j=ny1,ny2
1283	do i=nx1,nx2
1284
1285
1286	if (ghost_x(i,j)) then ! noeuds fantomes
1287	Tu(i,j,:,:)=0.
1288	Tv(i,j,:,:)=0.
1289	Tu(i,j,0,0)=1.
1290	opposx(i,j)=-3.33333
1291	end if
1292
1293	if (ghost_y(i,j)) then ! noeuds fantomes
1294	Su(i,j,:,:)=0.
1295	Sv(i,j,:,:)=0.
1296	Sv(i,j,0,0)=1.
1297	opposy(i,j)=-3.33333
1298	end if
1299
1300	end do
1301	end do
1302
1303	return
1304	end subroutine ghost_identite
1305	!------------------------------------------------------------------------------
1306
1307	subroutine ghost_remove ! met a 0 les Tuij ... qui corresponde a des ghost
1308
1309	do j=ny1,ny2
1310	do i=nx1,nx2
1311
1312	!u_ghost: if (ok_umat(i,j)) then
1313
1314
1315	ilmin=max(-1,nx1+1-i) ! pour avoir i+il entre nx1+1 et nx2
1316	ilmax=min(1,nx2-i)
1317
1318	jlmin=max(-1,ny1-j) ! pour avoir j+jl entre ny1 et ny2
1319	jlmax=min(1,ny2-j)
1320
1321
1322
1323	do jl = jlmin , jlmax ! balaye tous les coefficients d'un noeud U de -1 a 1
1324	do il = ilmin , ilmax ! de -1 a 1
1325
1326	if (ghost_x(i+il,j+jl)) then
1327	Tu(i,j,il,jl)=0.
1328	Su(i,j,il,jl)=0.
1329	end if
1330
1331	end do
1332	end do
1333	! end if u_ghost
1334
1335	!v_ghost: if (ok_vmat(i,j)) then
1336
1337	ilmin=max(-1,nx1-i) ! pour avoir i+il entre nx1 et nx2
1338	ilmax=min(1,nx2-i)
1339
1340	jlmin=max(-1,ny1+1-j) ! pour avoir j+jl entre ny1+1 et ny2
1341	jlmax=min(1,ny2-j)
1342
1343	do jl = jlmin , jlmax ! balaye tous les coefficients d'un noeud U de -1 a 1
1344	do il = ilmin , ilmax ! de -1 a 1
1345	if (ghost_y(i+il,j+jl)) then
1346	Tv(i,j,il,jl)=0.
1347	Sv(i,j,il,jl)=0.
1348	end if
1349	end do
1350	end do
1351
1352	! end if v_ghost
1353
1354	end do
1355	end do
1356	end subroutine ghost_remove
1357	!----------------------------------------------------------------------------------------
1358	subroutine calc_beta(uxgiven,uygiven) ! calcule betamx et betamy
1359	! a partir du champ de vitesse
1360	! uxprec, uyprec
1361	implicit none
1362	real, dimension(nx,ny) :: uxgiven ! prescribed velocity
1363	real, dimension(nx,ny) :: uygiven !
1364	real :: maxbeta = 1.e6 ! beta will take maxbeta value when
1365	real :: limgliss = 0.1 ! ugiven < limgliss (maxbeta in Pa an /m)
1366	real :: betamin = 10. ! minimum value on grounded point (in Pa an /m)
1367
1368	debug_3D(:,:,69)=uxgiven(:,:)
1369	debug_3D(:,:,70)=uygiven(:,:)
1370
1371
1372	if (itracebug.eq.1) call tracebug(' Subroutine calc_beta')
1373
1374	do j=ny1+1,ny2-1
1375	do i=nx1+1,nx2-1
1376
1377	x_calc: if (flgzmx(i,j)) then
1378
1379	betamx(i,j)=-opposx(i,j)
1380	debug_3D(i,j,35) = opposx(i,j)/dx/dx ! driving stress ou pression hydrostatique (Pa)
1381	do jl=-1,1
1382	do il=-1,1
1383	betamx(i,j) = betamx(i,j)+(Tu(i,j,il,jl)*uxgiven(i+il,j+jl) &
1384	+Tv(i,j,il,jl)*uygiven(i+il,j+jl))
1385
1386	debug_3D(i,j,82)=debug_3D(i,j,82) + (Tu(i,j,il,jl)*uxgiven(i+il,j+jl) &
1387	+Tv(i,j,il,jl)*uygiven(i+il,j+jl))
1388
1389	end do
1390	end do
1391
1392	debug_3D(i,j,82)=debug_3D(i,j,82)/dx/dx ! longitudinal stress computed from velocities
1393
1394
1395	! fonctionne même quand frotmx n'est pas nul
1396
1397	betamx(i,j) = betamx(i,j) + frotmx(i,j)dx2 uxgiven(i,j)
1398
1399	if (abs(uxgiven(i,j)).gt.limgliss) then
1400	betamx(i,j)=betamx(i,j)/dx2/uxgiven(i,j)
1401	betamx(i,j)=betamx(i,j)*drag_mx(i,j)
1402
1403	else
1404	if (flotmx(i,j)) then
1405	betamx(i,j)=0.
1406	else
1407	betamx(i,j)= maxbeta ! il faudrait mettre la limitation a l'aide du driving stress
1408	end if
1409	endif
1410
1411	else ! les noeuds non flgzmx sont d'office a maxbeta
1412	betamx(i,j)=maxbeta
1413	end if x_calc
1414
1415	y_calc: if (flgzmy(i,j)) then
1416
1417	betamy(i,j)=-opposy(i,j)
1418	debug_3D(i,j,36) = opposy(i,j)/dx/dx ! driving stress ou pression hydrostatique (Pa)
1419	do jl=-1,1
1420	do il=-1,1
1421
1422	betamy(i,j) = betamy(i,j)+(Su(i,j,il,jl)*uxgiven(i+il,j+jl) &
1423	+Sv(i,j,il,jl)*uygiven(i+il,j+jl))
1424
1425	debug_3D(i,j,83)=debug_3D(i,j,83) + (Su(i,j,il,jl)*uxgiven(i+il,j+jl) &
1426	+Sv(i,j,il,jl)*uygiven(i+il,j+jl))
1427
1428
1429	end do
1430	end do
1431
1432	debug_3D(i,j,83)=debug_3D(i,j,83) /dx/dx ! longitudinal stress computed from velocities
1433
1434	! fonctionne même quand frotmx n'est pas nul
1435
1436	betamy(i,j) = betamy(i,j) + frotmy(i,j)dx2 uygiven(i,j)
1437
1438	if (abs(uygiven(i,j)).gt.limgliss) then
1439	betamy(i,j)=betamy(i,j)/dx2/uygiven(i,j)
1440	betamx(i,j)=betamx(i,j)*drag_my(i,j)
1441	else
1442	if (flotmy(i,j)) then
1443	betamy(i,j)=0.
1444	else
1445	betamy(i,j)= maxbeta ! il faudrait mettre la limitation a l'aide du driving stress
1446	end if
1447	endif
1448
1449	else
1450	betamy(i,j)=maxbeta ! les noeuds non flgzmx sont d'office a maxbeta
1451	end if y_calc
1452
1453	end do
1454	end do
1455
1456	! loop to spread negative beta on the neighbours
1457
1458	do j=2,ny-1
1459	do i=2,nx-1
1460	if (betamx(i,j).lt.0.) then
1461	betamx(i-1,j) = max(betamin, (betamx(i-1,j) - betamx(i,j)) * 0.5)
1462	betamx(i+1,j) = max(betamin,(betamx(i+1,j) - betamx(i,j)) * 0.5)
1463	betamx(i,j) = max(betamin,betamx(i,j))
1464	end if
1465	if (betamy(i,j).lt.0.) then
1466	betamy(i,j-1) = max(betamin,betamy(i,j-1) - betamy(i,j) * 0.5)
1467	betamy(i,j+1) = max(betamin,betamy(i,j+1) - betamy(i,j) * 0.5)
1468	betamy(i,j) = max(betamin,betamy(i,j))
1469	end if
1470	end do
1471	end do
1472
1473	! average
1474	beta_centre(:,:)=0.
1475	do j=2,ny-1
1476	do i=2,nx-1
1477	beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
1478	+ (betamy(i,j)+betamy(i,j+1)))*0.25
1479	end do
1480	end do
1481
1482
1483
1484	! les noeuds negatifs veulent dire qu'il faudrait ajouter un moteur pour aller aussi vite
1485	! que les vitesses de bilan.
1486	! il faudrait repartir plutot que mettre a 0
1487
1488	! a mettre ailleurs
1489	! betamx(:,:)=max(betamx(:,:),0.)
1490	! betamy(:,:)=max(betamy(:,:),0.)
1491	! beta_centre(:,:)=max(beta_centre(:,:),0.)
1492
1493	! iter_beta=2
1494
1495
1496	if (itracebug.eq.1) call tracebug(' Fin calc_beta')
1497	end subroutine calc_beta
1498
1499	end subroutine rempli_L2

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

Download in other formats: