source: branches/GRISLIv3/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90 @ 467

Last change on this file since 467 was 467, checked in by aquiquet, 4 months ago

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