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

Last change on this file since 180 was 180, checked in by aquiquet, 6 years ago

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