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

Last change on this file since 4 was 4, checked in by dumas, 10 years ago

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