Index: /trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90
===================================================================
--- /trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90	(revision 73)
+++ /trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90	(revision 74)
@@ -41,4 +41,5 @@
 !-----------------------------------------------------------------------
 
+!$ USE OMP_LIB
 use module3d_phy
 use remplimat_declar        ! les tableaux tuij, ....
@@ -104,4 +105,6 @@
 
 !-----------------------------
+!$OMP PARALLEL
+!$OMP WORKSHARE
 Tu(:,:,:,:) = 0. ; Tv(:,:,:,:) = 0. ; Su(:,:,:,:) = 0. ; Sv(:,:,:,:) = 0.
 opposx(:,:) = 0. ; opposy(:,:) = 0.
@@ -117,5 +120,6 @@
 Hmx_oppos(:,:) = Hmx(:,:)
 Hmy_oppos(:,:) = Hmy(:,:)
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 ! calcul de Hmx_oppos et Hmy_oppos dans le cas de fronts
@@ -164,6 +168,10 @@
  call limit_frotm
 else
+   !$OMP PARALLEL
+   !$OMP WORKSHARE
    frotmx(:,:)= betamx(:,:)
    frotmy(:,:)= betamy(:,:)
+   !$OMP END WORKSHARE
+   !$OMP END PARALLEL
 end if
 
@@ -272,5 +280,6 @@
 ! pour avoir la diagonale=1 a faire imperativement APRES okmat
 !--------------------------------------------------------------
-
+!$OMP PARALLEL
+!$OMP DO PRIVATE(scal)
 do j=ny1,ny2
    do i=nx1,nx2
@@ -319,5 +328,6 @@
    end do
 end do
-
+!$OMP END DO
+!$OMP END PARALLEL
 
 ! mise a identite des noeuds fantomes
@@ -329,7 +339,6 @@
 
 
-
-debug_3d(:,:,35) = opposx(:,:)
-debug_3d(:,:,36) = opposy(:,:)
+!debug_3d(:,:,35) = opposx(:,:)
+!debug_3d(:,:,36) = opposy(:,:)
 
 
@@ -349,5 +358,6 @@
 
 ! remise a 0 des noeuds fantomes
-
+!$OMP PARALLEL
+!$OMP WORKSHARE
 where ((Hmx(nx1:nx2,ny1:ny2).le.1.).and.(flgzmx(nx1:nx2,ny1:ny2)))
    uxnew(nx1:nx2,ny1:ny2)=0.
@@ -357,5 +367,6 @@
    uynew(nx1:nx2,ny1:ny2)=0.
 end where
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 ! call graphique_L2(kl,ku,nx1,nx2,ny1,ny2,imx(nx1:nx2,ny1:ny2),imy(nx1:nx2,ny1:ny2))   
@@ -404,4 +415,6 @@
 
 subroutine limit_frotm
+
+!$USE OMP_LIB
 !-------------------------
 ! limite la valeur de frotmx a betamax
@@ -411,4 +424,6 @@
 if (itracebug.eq.1) write(Num_tracebug,*)'betamax = ',betamax
 
+!$OMP PARALLEL
+!$OMP WORKSHARE
 where (flgzmx(:,:))
    frotmx(:,:)=min(abs(betamx(:,:)),betamax_2d(:,:))
@@ -422,4 +437,6 @@
    frotmy(:,:)=0
 end where
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 end subroutine limit_frotm
@@ -433,4 +450,6 @@
 
 subroutine rempli_Tuij                                    ! appelle tous les cas
+
+!$ USE OMP_LIB
 
 ! imx(i,j)=0 -> ne pas traiter cette vitesse
@@ -464,5 +483,6 @@
 count_line=1                                 ! pour numeroter les lignes
 !------------------------------------------------------------------
-
+!$OMP PARALLEL
+!$OMP DO
 lignes_UV: do j=ny1,ny2 
    do i=nx1,nx2
@@ -470,4 +490,5 @@
       if (i.gt.nx1) then                      ! Vitesses U
          ligu_L2(i,j)=count_line
+         !$OMP ATOMIC
          count_line=count_line+1
           
@@ -480,32 +501,33 @@
 
          case(1)                              ! vitesse imposee
-            call vel_U_presc     
+            call vel_U_presc(uxprec(i,j),Tu(i,j,0,0),opposx(i,j))    
 
          case(2)                              ! cas general
-            call vel_U_general
+            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), &
+            Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
 
          case(-1)                             ! bord sud
-            call vel_U_South          
+            call vel_U_South(Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))          
 
          case(-2)                             ! bord Est
-            call vel_U_East
+            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))
        
          case(-3)                             ! bord nord
-            call vel_U_North   
+            call vel_U_North(Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))   
 
          case(-4)                             ! bord West
-            call vel_U_West
+            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))
 
          case(-12)                            ! coin SE
-            call vel_U_SE
+            call vel_U_SE(Tu(i,j,:,:),opposx(i,j))
 
          case(-23)                            ! coin NE
-            call vel_U_NE
+            call vel_U_NE(Tu(i,j,:,:),opposx(i,j))
 
          case(-34)                            ! coin NW
-            call vel_U_NW
+            call vel_U_NW(Tu(i,j,:,:),opposx(i,j))
 
          case(-41)                            ! coin SW
-            call vel_U_SW
+            call vel_U_SW(Tu(i,j,:,:),opposx(i,j))
 
          end select case_imx
@@ -515,4 +537,5 @@
 
          ligv_L2(i,j)=count_line
+         !$OMP ATOMIC
          count_line=count_line+1
 
@@ -525,33 +548,34 @@
 
          case(1)                              ! vitesse imposee
-            call vel_V_presc
+            call vel_V_presc(uyprec(i,j),Sv(i,j,0,0),opposy(i,j))
             
          case(2)                              ! cas general
-            call vel_V_general
+            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), &
+            Su(i,j,:,:),Sv(i,j,:,:),opposy(i,j))
             
          case(-1)                             ! bord sud
-            call vel_V_South             
+            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))
             
          case(-2)                             ! bord Est
-            call vel_V_East
+            call vel_V_East(Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
             
          case(-3)                             ! bord nord
-            call vel_V_North          
+            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))
             
          case(-4)                             ! bord West
-            call vel_V_West
+            call vel_V_West(Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
             
          case(-12)                            ! coin SE
-            call vel_V_SE
+            call vel_V_SE(Sv(i,j,:,:),opposy(i,j)) 
             
          case(-23)                            ! coin NE
-            call vel_V_NE
+            call vel_V_NE(Sv(i,j,:,:),opposy(i,j))
             
          case(-34)                            ! coin NW
-            call vel_V_NW
+            call vel_V_NW(Sv(i,j,:,:),opposy(i,j)) 
             
          case(-41)                            ! coin SW
 
-            call vel_V_SW
+            call vel_V_SW(Sv(i,j,:,:),opposy(i,j))
             
          end select case_imy
@@ -559,4 +583,6 @@
    end do
 end do lignes_UV
+!$OMP END DO
+!$OMP END PARALLEL
 
 end subroutine rempli_Tuij
@@ -570,38 +596,56 @@
 !------------------------------------------------------------------
 
-subroutine vel_U_presc     !                     ! vitesse imposee
-
-  Tu(i,j,0,0)=1.
-  opposx(i,j)=uxprec(i,j)
+subroutine vel_U_presc(uxprec,Tu,opposx)     !                     ! vitesse imposee
+
+implicit none
+
+real,intent(in)  :: uxprec
+real,intent(out) :: Tu
+real,intent(out) :: opposx
+
+  Tu=1.
+  opposx=uxprec
 
 end subroutine vel_U_presc
 !------------------------------------------------------------------
                    
-subroutine vel_U_general                   ! cas general
-
-
-   beta=frotmx(i,j)*dx2                       !  Terme en u(i,j)
-
-
-
-   Tu(i,j,0,0)  = -4.*pvi(i,j)-4.*pvi(i-1,j)- pvm(i,j+1)-pvm(i,j)-beta
-
-   Tu(i,j,-1,0) =  4.*pvi(i-1,j)              !  Terme en u(i-1,j)
-
-   Tu(i,j,1,0)  =  4.*pvi(i,j)                !  Terme en u(i+1,j)
-
-   Tu(i,j,0,-1) =  pvm(i,j)                   !  Terme en u(i,j-1)
-
-   Tu(i,j,0,1)  =  pvm(i,j+1)                 !  Terme en u(i,j+1)
-
-   Tv(i,j,0,0)  = -2.*pvi(i,j)-pvm(i,j)       !  Terme en v(i,j)
-
-   Tv(i,j,-1,0) =  2.*pvi(i-1,j)+pvm(i,j)     !  Terme en v(i-1,j)
-
-   Tv(i,j,-1,1) = -2.*pvi(i-1,j)-pvm(i,j+1)   !  Terme en v(i-1,j+1)
-
-   Tv(i,j,0,1)  =  2.*pvi(i,j)+pvm(i,j+1)     !  Terme en v(i,j+1)
-
-   opposx(i,j)  =  rog*hmx_oppos(i,j)*sdx(i,j)*dx2  !  Terme en opposx(i,j)
+subroutine vel_U_general(frotmx,dx2,pvi,pvi_imoinsun,pvm,pvm_jplusun,sdx,rog,hmx_oppos,Tu,Tv,opposx)                   ! cas general
+
+implicit none
+
+real,intent(in) :: frotmx
+real,intent(in) :: dx2
+real,intent(in) :: pvi
+real,intent(in) :: pvi_imoinsun
+real,intent(in) :: pvm
+real,intent(in) :: pvm_jplusun
+real,intent(in) :: sdx
+real,intent(in) :: rog
+real,intent(in) :: hmx_oppos
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,dimension(-1:1,-1:1),intent(inout) :: Tv
+real,intent(out) :: opposx
+
+   beta=frotmx*dx2                       !  Terme en u(i,j)
+
+   Tu(0,0)  = -4.*pvi - 4.*pvi_imoinsun - pvm_jplusun - pvm - beta
+
+   Tu(-1,0) =  4.*pvi_imoinsun              !  Terme en u(i-1,j)
+
+   Tu(1,0)  =  4.*pvi                !  Terme en u(i+1,j)
+
+   Tu(0,-1) =  pvm                   !  Terme en u(i,j-1)
+
+   Tu(0,1)  =  pvm_jplusun                 !  Terme en u(i,j+1)
+
+   Tv(0,0)  = -2.*pvi - pvm       !  Terme en v(i,j)
+
+   Tv(-1,0) =  2.*pvi_imoinsun + pvm     !  Terme en v(i-1,j)
+
+   Tv(-1,1) = -2.*pvi_imoinsun - pvm_jplusun   !  Terme en v(i-1,j+1)
+
+   Tv(0,1)  =  2.*pvi + pvm_jplusun     !  Terme en v(i,j+1)
+
+   opposx   =  rog*hmx_oppos*sdx*dx2  !  Terme en opposx(i,j)
 
    return
@@ -611,11 +655,16 @@
 ! voir explications dans vel_U_West
 
-subroutine vel_U_South                       ! bord sud non cisaillement           
-
-  Tu(i,j,0,0)  =  1.
-  Tu(i,j,0,1)  = -1.
-  Tv(i,j,0,1)  = -1.
-  Tv(i,j,-1,1) =  1.
-  opposx(i,j)  =  0.
+subroutine vel_U_South(Tu,Tv,opposx)       ! bord sud non cisaillement           
+
+implicit none                       
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,dimension(-1:1,-1:1),intent(inout) :: Tv
+real, intent(out) :: opposx
+
+  Tu(0,0)  =  1.
+  Tu(0,1)  = -1.
+  Tv(0,1)  = -1.
+  Tv(-1,1) =  1.
+  opposx   =  0.
 
    return
@@ -623,11 +672,16 @@
 !------------------------------------------------------------------
 
-subroutine vel_U_North                       ! bord nord   non cisaillement       
-
-  Tu(i,j,0,0)  =  1.
-  Tu(i,j,0,-1) = -1.
-  Tv(i,j,0,0)  =  1.
-  Tv(i,j,-1,0) = -1.
-  opposx(i,j)  =  0.
+subroutine vel_U_North(Tu,Tv,opposx)     ! bord nord   non cisaillement       
+
+implicit none                       
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,dimension(-1:1,-1:1),intent(inout) :: Tv
+real, intent(out) :: opposx
+
+  Tu(0,0)  =  1.
+  Tu(0,-1) = -1.
+  Tv(0,0)  =  1.
+  Tv(-1,0) = -1.
+  opposx   =  0.
 
    return
@@ -636,11 +690,24 @@
 ! voir explications dans vel_U_West
 
-subroutine vel_U_East                       ! bord Est   pression
+subroutine vel_U_East(pvi,pvi_imoinsun,rog,H_imoinsun,rowg,slv_imoinsun,B_imoinsun,dx,Tu,Tv,opposx) ! bord Est   pression
+
+implicit none
+real,intent(in) :: pvi
+real,intent(in) :: pvi_imoinsun
+real,intent(in) :: rog
+real,intent(in) :: H_imoinsun
+real,intent(in) :: rowg
+real,intent(in) :: slv_imoinsun
+real,intent(in) :: B_imoinsun
+real,intent(in) :: dx                       
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,dimension(-1:1,-1:1),intent(inout) :: Tv
+real, intent(out) :: opposx
        
-  Tu(i,j,0,0)   =  4.*pvi(i-1,j)
-  Tu(i,j,-1,0)  = -4.*pvi(i-1,j)
-  Tv(i,j,0,1)   =  2.*pvi(i,j)
-  Tv(i,j,0,0)   = -2.*pvi(i,j)
-  opposx(i,j)   =  0.5*(rog*H(i-1,j)**2-rowg*(max(slv(i-1,j)-B(i-1,j),0.))**2)*dx
+  Tu(0,0)   =  4.*pvi_imoinsun
+  Tu(-1,0)  = -4.*pvi_imoinsun
+  Tv(0,1)   =  2.*pvi
+  Tv(0,0)   = -2.*pvi
+  opposx    =  0.5*(rog*H_imoinsun**2-rowg*(max(slv_imoinsun-B_imoinsun,0.))**2)*dx
 
    return
@@ -648,12 +715,25 @@
 !------------------------------------------------------------------
 
-subroutine vel_U_West                        ! bord West pression
+subroutine vel_U_West(pvi,pvi_imoinsun,rog,H,rowg,slv,B,dx,Tu,Tv,opposx) ! bord West pression
                                              ! tous les coef * -1
+implicit none
+real,intent(in) :: pvi
+real,intent(in) :: pvi_imoinsun
+real,intent(in) :: rog
+real,intent(in) :: H
+real,intent(in) :: rowg
+real,intent(in) :: slv
+real,intent(in) :: B
+real,intent(in) :: dx                       
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,dimension(-1:1,-1:1),intent(inout) :: Tv
+real, intent(out) :: opposx
+
 
 !  Tu(i,j,0,0)    =  4.*pvi(i-1,j)
 !  Tu(i,j,1,0)    = -4.*pvi(i-1,j)
 
-  Tu(i,j,0,0)    =  4.*pvi(i,j)              ! quand le bord est epais, il faut prendre 
-  Tu(i,j,1,0)    = -4.*pvi(i,j)              ! le noeud du shelf pas celui en H=1
+  Tu(0,0)    =  4.*pvi              ! quand le bord est epais, il faut prendre 
+  Tu(1,0)    = -4.*pvi              ! le noeud du shelf pas celui en H=1
                                              ! s'il est fin, pas de difference
 
@@ -661,7 +741,7 @@
 !!$  Tv(i-1,j,-1,1) = -2.*pvi(i-1,j)
 
-  Tv(i,j,-1,0) =  2.*pvi(i-1,j)
-  Tv(i,j,-1,1) = -2.*pvi(i-1,j)
- opposx(i,j)    = -0.5*(rog*H(i,j)**2-rowg*(max(slv(i,j)-B(i,j),0.))**2)*dx
+  Tv(-1,0) =  2.*pvi_imoinsun
+  Tv(-1,1) = -2.*pvi_imoinsun
+ opposx    = -0.5*(rog*H**2-rowg*(max(slv-B,0.))**2)*dx
 
  
@@ -674,8 +754,13 @@
 !------------------------------------------------------------------
 
-subroutine vel_V_presc                     ! vitesse imposee
-
-  Sv(i,j,0,0)=1.
-  opposy(i,j)=uyprec(i,j)
+subroutine vel_V_presc(uyprec,Sv,opposy)          ! vitesse imposee
+
+implicit none
+real,intent(in)  :: uyprec
+real,intent(out) :: Sv
+real,intent(out) :: opposy
+
+  Sv=1.
+  opposy=uyprec
 
    return
@@ -683,27 +768,43 @@
 !------------------------------------------------------------------
                     
-subroutine vel_V_general                      ! cas general
-
-
-   beta          =  frotmy(i,j)*dx2           !  Terme en v(i,j)
-   Sv(i,j,0,0)   = -4.*pvi(i,j)-4.*pvi(i,j-1)- pvm(i+1,j)-pvm(i,j)-beta
-
-   Sv(i,j,0,-1)  =  4.*pvi(i,j-1)             !  Terme en v(i,j-1)
-
-   Sv(i,j,0,1)   =  4.*pvi(i,j)               !  Terme en v(i,j+1)
-
-   Sv(i,j,-1,0)  =  pvm(i,j)                  !  Terme en v(i-1,j)
-
-   Sv(i,j,1,0)   =  pvm(i+1,j)                !  Terme en v(i+1,j)
-
-   Su(i,j,0,0)   = -2.*pvi(i,j)-pvm(i,j)      !  Terme en u(i,j)
-
-   Su(i,j,0,-1)  =  2.*pvi(i,j-1)+pvm(i,j)    !  Terme en u(i,j-1)
-
-   Su(i,j,1,-1)  = -2.*pvi(i,j-1)-pvm(i+1,j)  !  Terme en u(i+1,j-1)
-
-   Su(i,j,1,0)   =  2.*pvi(i,j)+pvm(i+1,j)    !  Terme en u(i+1,j)
-
-   opposy(i,j)   =  rog*hmy_oppos(i,j)*sdy(i,j)*dx2 !  Terme en opposy(i,j)
+subroutine vel_V_general(frotmy,dx2,pvi,pvi_jmoinsun,pvm,pvm_iplusun,sdy,rog,hmy_oppos,Su,Sv,opposy)! cas general
+
+implicit none
+real,intent(in) :: frotmy
+real,intent(in) :: dx2
+real,intent(in) :: pvi
+real,intent(in) :: pvi_jmoinsun
+real,intent(in) :: pvm
+real,intent(in) :: pvm_iplusun
+real,intent(in) :: sdy
+real,intent(in) :: rog
+real,intent(in) :: hmy_oppos
+real,dimension(-1:1,-1:1),intent(inout) :: Su
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,intent(out) :: opposy
+
+real :: beta
+
+
+   beta          =  frotmy*dx2           !  Terme en v(i,j)
+   Sv(0,0)   = -4.*pvi - 4.*pvi_jmoinsun - pvm_iplusun - pvm-beta
+
+   Sv(0,-1)  =  4.*pvi_jmoinsun             !  Terme en v(i,j-1)
+
+   Sv(0,1)   =  4.*pvi               !  Terme en v(i,j+1)
+
+   Sv(-1,0)  =  pvm                  !  Terme en v(i-1,j)
+
+   Sv(1,0)   =  pvm_iplusun                !  Terme en v(i+1,j)
+
+   Su(0,0)   = -2.*pvi - pvm      !  Terme en u(i,j)
+
+   Su(0,-1)  =  2.*pvi_jmoinsun + pvm    !  Terme en u(i,j-1)
+
+   Su(1,-1)  = -2.*pvi_jmoinsun - pvm_iplusun  !  Terme en u(i+1,j-1)
+
+   Su(1,0)   =  2.*pvi + pvm_iplusun    !  Terme en u(i+1,j)
+
+   opposy    =  rog*hmy_oppos*sdy*dx2 !  Terme en opposy(i,j)
 
    return
@@ -712,11 +813,16 @@
 !------------------------------------------------------------------
 
-subroutine vel_V_West                       ! bord West non cisaillement           
-
-  Sv(i,j,0,0)  =  1.
-  Sv(i,j,1,0)  = -1.
-  Su(i,j,1,0)  = -1.
-  Su(i,j,1,-1) =  1.
-  opposy(i,j)  =  0.
+subroutine vel_V_West(Sv,Su,opposy)     ! bord West non cisaillement           
+
+implicit none
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,dimension(-1:1,-1:1),intent(inout) :: Su
+real,intent(out) :: opposy
+
+  Sv(0,0)  =  1.
+  Sv(1,0)  = -1.
+  Su(1,0)  = -1.
+  Su(1,-1) =  1.
+  opposy   =  0.
 
    return
@@ -725,11 +831,16 @@
 
 
-subroutine vel_V_East                       ! bord Est   non cisaillement       
-
-  Sv(i,j,0,0)  =  1.
-  Sv(i,j,-1,0) = -1.
-  Su(i,j,0,0)  =  1.
-  Su(i,j,0,-1) = -1.
-  opposy(i,j)  =  0.
+subroutine vel_V_East(Sv,Su,opposy)         ! bord Est   non cisaillement       
+
+implicit none
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,dimension(-1:1,-1:1),intent(inout) :: Su
+real,intent(out) :: opposy
+
+  Sv(0,0)  =  1.
+  Sv(-1,0) = -1.
+  Su(0,0)  =  1.
+  Su(0,-1) = -1.
+  opposy   =  0.
 
    return
@@ -739,11 +850,24 @@
 ! voir explications dans vel_U_West
 
-subroutine vel_V_North                       ! bord Nord   pression
+subroutine vel_V_North(pvi,pvi_jmoinsun,rog,H_imoinsun,rowg,slv_imoinsun,B_imoinsun,dx,Sv,Su,opposy) ! bord Nord   pression
+
+implicit none
+real,intent(in) :: pvi
+real,intent(in) :: pvi_jmoinsun
+real,intent(in) :: rog
+real,intent(in) :: H_imoinsun
+real,intent(in) :: rowg
+real,intent(in) :: slv_imoinsun
+real,intent(in) :: B_imoinsun
+real,intent(in) :: dx
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,dimension(-1:1,-1:1),intent(inout) :: Su
+real,intent(out) :: opposy
        
-  Sv(i,j,0,0)  =  4.*pvi(i,j-1)
-  Sv(i,j,0,-1) = -4.*pvi(i,j-1)
-  Su(i,j,1,0)  =  2.*pvi(i,j)
-  Su(i,j,0,0)  = -2.*pvi(i,j)
-  opposy(i,j)  = 0.5*(rog*H(i-1,j)**2-rowg*(max(slv(i-1,j)-B(i-1,j),0.))**2)*dx
+  Sv(0,0)  =  4.*pvi_jmoinsun
+  Sv(0,-1) = -4.*pvi_jmoinsun
+  Su(1,0)  =  2.*pvi
+  Su(0,0)  = -2.*pvi
+  opposy   = 0.5*(rog*H_imoinsun**2-rowg*(max(slv_imoinsun-B_imoinsun,0.))**2)*dx
 
    return
@@ -752,12 +876,24 @@
 ! voir explications dans vel_U_West
 
-subroutine vel_V_South                        ! bord sud  pression
+subroutine vel_V_South(pvi,pvi_jmoinsun,rog,H,rowg,slv,B,dx,Sv,Su,opposy) ! bord sud  pression
                                               ! tous les coef * -1
-
-  Sv(i,j,0,0)  = 4.*pvi(i,j)
-  Sv(i,j,0,1)  =-4.*pvi(i,j)
-  Su(i,j,0,-1) = 2.*pvi(i,j-1)
-  Su(i,j,1,-1) =-2.*pvi(i,j-1)
-  opposy(i,j)  = -0.5*(rog*H(i,j)**2-rowg*(max(slv(i,j)-B(i,j),0.))**2)*dx
+implicit none
+real,intent(in) :: pvi
+real,intent(in) :: pvi_jmoinsun
+real,intent(in) :: rog
+real,intent(in) :: H
+real,intent(in) :: rowg
+real,intent(in) :: slv
+real,intent(in) :: B
+real,intent(in) :: dx
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,dimension(-1:1,-1:1),intent(inout) :: Su
+real,intent(out) :: opposy
+
+  Sv(0,0)  = 4.*pvi
+  Sv(0,1)  =-4.*pvi
+  Su(0,-1) = 2.*pvi_jmoinsun
+  Su(1,-1) =-2.*pvi_jmoinsun
+  opposy   = -0.5*(rog*H**2-rowg*(max(slv-B,0.))**2)*dx
 
    return
@@ -780,9 +916,13 @@
 !Coin SE
 !----------
-subroutine vel_U_SE                        
-
-  Tu(i,j,0,0) =  1.    !*pvi(i,j)
-  Tu(i,j,0,1) = -1.    !*pvi(i,j)
-  opposx(i,j) =  0.
+subroutine vel_U_SE(Tu,opposx)                        
+
+implicit none
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,intent(out) :: opposx
+
+  Tu(0,0) =  1.    !*pvi(i,j)
+  Tu(0,1) = -1.    !*pvi(i,j)
+  opposx  =  0.
 
    return
@@ -790,9 +930,13 @@
 !------------------------------------------------------------------
 
-subroutine vel_V_SE                        
-
-  Sv(i,j,0,0) =  1.   !*pvi(i,j)
-  Sv(i,j,-1,0)= -1.   !*pvi(i,j)
-  opposy(i,j) =  0.
+subroutine vel_V_SE(Sv,opposy)                        
+
+implicit none
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,intent(out) :: opposy
+
+  Sv(0,0) =  1.   !*pvi(i,j)
+  Sv(-1,0)= -1.   !*pvi(i,j)
+  opposy  =  0.
 
    return
@@ -803,9 +947,13 @@
 !Coin SW
 !----------
-subroutine vel_U_SW                        
-
-  Tu(i,j,0,0) =  1.   ! *pvi(i,j)
-  Tu(i,j,0,1) = -1.   ! *pvi(i,j)
-  opposx(i,j) =  0.
+subroutine vel_U_SW(Tu,opposx)                        
+
+implicit none
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,intent(out) :: opposx
+
+  Tu(0,0) =  1.   ! *pvi(i,j)
+  Tu(0,1) = -1.   ! *pvi(i,j)
+  opposx  =  0.
 
    return
@@ -813,9 +961,13 @@
 !------------------------------------------------------------------
 
-subroutine vel_V_SW                        
-
-  Sv(i,j,0,0) =  1.   !*pvi(i,j)
-  Sv(i,j,1,0) = -1.   !*pvi(i,j)
-  opposy(i,j) =  0.
+subroutine vel_V_SW(Sv,opposy)
+
+implicit none
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,intent(out) :: opposy                    
+
+  Sv(0,0) =  1.   !*pvi(i,j)
+  Sv(1,0) = -1.   !*pvi(i,j)
+  opposy  =  0.
 
    return
@@ -825,9 +977,13 @@
 !Coin NE
 !----------
-subroutine vel_U_NE                        
-
-  Tu(i,j,0,0)  =  1.  ! *pvi(i,j)
-  Tu(i,j,0,-1) = -1.  ! *pvi(i,j)
-  opposx(i,j)  =  0.
+subroutine vel_U_NE(Tu,opposx)
+
+implicit none
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,intent(out) :: opposx                     
+
+  Tu(0,0)  =  1.  ! *pvi(i,j)
+  Tu(0,-1) = -1.  ! *pvi(i,j)
+  opposx   =  0.
 
    return
@@ -835,9 +991,13 @@
 !------------------------------------------------------------------
 
-subroutine vel_V_NE                        
-
-  Sv(i,j,0,0)  =  1.  ! *pvi(i,j)
-  Sv(i,j,-1,0) = -1.  ! *pvi(i,j)
-  opposy(i,j)  =  0.
+subroutine vel_V_NE(Sv,opposy)                        
+
+implicit none
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,intent(out) :: opposy
+
+  Sv(0,0)  =  1.  ! *pvi(i,j)
+  Sv(-1,0) = -1.  ! *pvi(i,j)
+  opposy   =  0.
 
    return
@@ -847,9 +1007,13 @@
 !Coin NW
 !----------
-subroutine vel_U_NW                        
-
-  Tu(i,j,0,0)  =  1.  ! *pvi(i,j)
-  Tu(i,j,0,-1) = -1.  ! *pvi(i,j)
-  opposx(i,j)  =  0.
+subroutine vel_U_NW(Tu,opposx)
+
+implicit none                        
+real,dimension(-1:1,-1:1),intent(inout) :: Tu
+real,intent(out) :: opposx 
+
+  Tu(0,0)  =  1.  ! *pvi(i,j)
+  Tu(0,-1) = -1.  ! *pvi(i,j)
+  opposx   =  0.
 
    return
@@ -857,9 +1021,13 @@
 !------------------------------------------------------------------
 
-subroutine vel_V_NW                        
-
-  Sv(i,j,0,0)  =  1.  !*pvi(i,j)
-  Sv(i,j,1,0)  = -1.  !*pvi(i,j)
-  opposy(i,j)  =  0.
+subroutine vel_V_NW(Sv,opposy)                        
+
+implicit none                        
+real,dimension(-1:1,-1:1),intent(inout) :: Sv
+real,intent(out) :: opposy 
+
+  Sv(0,0)  =  1.  !*pvi(i,j)
+  Sv(1,0)  = -1.  !*pvi(i,j)
+  opposy   =  0.
 
    return
Index: /trunk/SOURCES/litho-0.4.f90
===================================================================
--- /trunk/SOURCES/litho-0.4.f90	(revision 73)
+++ /trunk/SOURCES/litho-0.4.f90	(revision 74)
@@ -37,63 +37,66 @@
 
  
-      subroutine litho
+subroutine litho
+  !$ USE OMP_LIB
+  USE module3D_phy
+  USE param_phy_mod
+  USE ISO_DECLAR ! module de declaration des variables specifiques a l'isostasie
 
-      USE module3D_phy
-      USE param_phy_mod
-      USE ISO_DECLAR ! module de declaration des variables specifiques a l'isostasie
+  implicit none
 
-      implicit none
-
-      INTEGER :: II,SOM1,SOM2
-      REAL, dimension(:,:), allocatable :: WLOC
+  INTEGER :: II,SOM1,SOM2
+  REAL, dimension(:,:), allocatable :: WLOC
 
 
-!----- allocation de WLOC  et de croix -----------
+  !----- allocation de WLOC  et de croix -----------
 
-      if (.not.allocated(WLOC)) THEN
-         allocate(WLOC(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err)
-         if (err/=0) then
-            print *,"Erreur a l'allocation du tableau WLOC",err
-            stop 4
-         end if
-      end if
+  if (.not.allocated(WLOC)) THEN
+     allocate(WLOC(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err)
+     if (err/=0) then
+	print *,"Erreur a l'allocation du tableau WLOC",err
+	stop 4
+     end if
+  end if
 
-!----- fin de l'allocation --------------
+  !----- fin de l'allocation --------------
 
-!     calcul de la deflexion par sommation des voisins appartenant
-!     au bloc de taille 2*LBLOC
-       som1=0.
-       som2=0.
+  !     calcul de la deflexion par sommation des voisins appartenant
+  !     au bloc de taille 2*LBLOC
+  som1=0.
+  som2=0.
 
-!      On somme aussi les contributions des points exterieurs au domaine
-!      lorsque la charge est due a l'ocean. On suppose alors  que
-!      ces points ont la meme charge que les limites
-
-boucleij : do J=1,NY
-           do I=1,NX
+  !      On somme aussi les contributions des points exterieurs au domaine
+  !      lorsque la charge est due a l'ocean. On suppose alors  que
+  !      ces points ont la meme charge que les limites
+  !$OMP PARALLEL
+  !$OMP DO PRIVATE(Wloc)
+  boucleij : do J=1,NY
+     do I=1,NX
 
 
-              W1(I,J)=0. 
-              ii=0
+	W1(I,J)=0. 
+ !   ii=0
 
-! Wloc : enfoncement centre sur ij du a chaque charges autour
-      Wloc(:,:) = We(:,:) * charge(i-lbloc:i+lbloc,j-lbloc:j+lbloc)
+ ! Wloc : enfoncement centre sur ij du a chaque charges autour
+	Wloc(:,:) = We(:,:) * charge(i-lbloc:i+lbloc,j-lbloc:j+lbloc)
 
 
-      W1(i,j) = sum (Wloc(:,:))   ! effet en ij de toutes les charges autour (distance < lbloc)
+	W1(i,j) = sum (Wloc(:,:))   ! effet en ij de toutes les charges autour (distance < lbloc)
 
-!    ------------------------ FIN DE L'INTEGRATION SUR LE PAVE LBLOC
-              som1=som1+W1(I,J)
-              som2=som2-charge(I,J)/ROMG
+ !    ------------------------ FIN DE L'INTEGRATION SUR LE PAVE LBLOC
+ !   som1=som1+W1(I,J)
+ !   som2=som2-charge(I,J)/ROMG
 
-      end do
-      end do  boucleij
-       
-!     write(6,*)'enfoncement total avec rigidite de la lithosphere:',
-!    &  som1
-!     write(6,*)'enfoncement total avec isostasie locale :',
-!    &  som2
+     end do
+  end do  boucleij
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  !     write(6,*)'enfoncement total avec rigidite de la lithosphere:',
+  !    &  som1
+  !     write(6,*)'enfoncement total avec isostasie locale :',
+  !    &  som2
 
 
-    end subroutine litho
-      
+end subroutine litho
+
Index: /trunk/SOURCES/taubed-0.3.f90
===================================================================
--- /trunk/SOURCES/taubed-0.3.f90	(revision 73)
+++ /trunk/SOURCES/taubed-0.3.f90	(revision 74)
@@ -43,4 +43,5 @@
 subroutine taubed()
 
+  !$USE OMP_LIB 
   USE module3D_phy
   USE param_phy_mod
@@ -52,6 +53,9 @@
   ! NLITH est defini dans isostasie et permet le choix du modele d'isostasie
 
+  
   if (NLITH.eq.1) then
      !       avec rigidite de la lithosphere
+     !$OMP PARALLEL
+     !$OMP DO
      do J=1,NY 
         do I=1,NX
@@ -65,28 +69,35 @@
         end do
      end do
-
+     !$OMP END DO
 
      ! il faut remplir CHARGE dans les parties a l'exterieur de la grille :
      ! a l'exterieur de la grille CHARGE est egale a la valeur sur le bord
-
+     !$OMP DO
      do J=1,NY
         CHARGE(1-LBLOC:0,J)=CHARGE(1,J)      ! bord W
         CHARGE(NX+1:NX+LBLOC,J)=CHARGE(NX,J) ! bord E
      end do
+     !$OMP END DO
+     !$OMP DO
      do I=1,NX
         CHARGE(I,1-LBLOC:0)=CHARGE(I,1)      ! bord S
         CHARGE(I,NY+1:NY+LBLOC)=CHARGE(I,NY) ! bord N
      end do
+     !$OMP END DO
 
      ! valeur dans les quatres coins exterieurs au domaine       
+	  !$OMP WORKSHARE
      CHARGE(1-LBLOC:0,1-LBLOC:0)=CHARGE(1,1)           ! coin SW
      CHARGE(1-LBLOC:0,NY+1:NY+LBLOC)=CHARGE(1,NY)      ! coin NW
      CHARGE(NX+1:NX+LBLOC,1-LBLOC:0)=CHARGE(NX,1)      ! coin SE
      CHARGE(NX+1:NX+LBLOC,NY+1:NY+LBLOC)=CHARGE(NX,NY) ! coin NE
-
+     !$OMP END WORKSHARE
+     !$OMP END PARALLEL 
      call litho
 
   else
      !     enfoncement local
+     !$OMP PARALLEL
+     !$OMP DO
      do J=1,NY
         do I=1,NX
@@ -100,7 +111,11 @@
         end do
      end do
+     !$OMP END DO
+     !$OMP END PARALLEL
   endif
 
   !     decroissance exponentielle de l'enfoncement
+  !$OMP PARALLEL
+  !$OMP DO
   do J=1,NY
      do I=1,NX
@@ -109,4 +124,6 @@
      end do
   end do
+  !$OMP END DO
+  !$OMP END PARALLEL
 
 end subroutine taubed
