Index: trunk/NEMO/OPA_SRC/SOL/sol_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/sol_oce.F90	(revision 311)
+++ trunk/NEMO/OPA_SRC/SOL/sol_oce.F90	(revision 312)
@@ -24,5 +24,5 @@
    !! ----------------------------------
    INTEGER , PUBLIC ::      & !!: namsol   elliptic solver / island / free surface
-      nsolv    =    1 ,     &  !: = 1/2/3 type of elliptic solver
+      nsolv    =    1 ,     &  !: = 1/2/3/4 type of elliptic solver
       nsol_arp =    0 ,     &  !: = 0/1 absolute/relative precision convergence test
       nmin     =  300 ,     &  !: minimum of iterations for the SOR solver
@@ -55,8 +55,8 @@
       rr               !: coefficient  =(rn,rn)
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,4) ::   &  !:
+   REAL(wp), PUBLIC, DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4) ::   &  !:
       gcp              !: barotropic matrix extra-diagonal elements
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
+   REAL(wp), PUBLIC, DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) ::   &  !:
       gcx, gcxb,    &  !: now, before solution of the elliptic equation
       gcdprc,       &  !: inverse diagonal preconditioning matrix
Index: trunk/NEMO/OPA_SRC/SOL/solmat.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solmat.F90	(revision 311)
+++ trunk/NEMO/OPA_SRC/SOL/solmat.F90	(revision 312)
@@ -17,4 +17,5 @@
    USE phycst          ! physical constants
    USE obc_oce         ! ocean open boundary conditions
+   USE lbclnk          ! lateral boudary conditions
    USE lib_mpp         ! distributed memory computing
    USE dynspg_rl
@@ -232,5 +233,15 @@
          gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:)
          gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:)
-         IF( nsolv == 2 )   gccd(:,:) = sor * gcp(:,:,2)
+         IF( ( nsolv == 2 ) .OR. ( nsolv == 2 ) )  gccd(:,:) = sor * gcp(:,:,2)
+
+         IF( nsolv == 4 ) THEN
+            CALL lbc_lnk_e( gcp   (:,:,1), c_solver_pt, 1. )   ! lateral boundary conditions
+            CALL lbc_lnk_e( gcp   (:,:,2), c_solver_pt, 1. )   ! lateral boundary conditions
+            CALL lbc_lnk_e( gcp   (:,:,3), c_solver_pt, 1. )   ! lateral boundary conditions
+            CALL lbc_lnk_e( gcp   (:,:,4), c_solver_pt, 1. )   ! lateral boundary conditions
+            CALL lbc_lnk_e( gcdprc(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions
+            CALL lbc_lnk_e( gcdmat(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions         
+            IF( npolj /= 0 ) CALL sol_exd( gcp , c_solver_pt ) ! switch northernelements
+         END IF
 
       ELSE
@@ -308,4 +319,104 @@
       
    END SUBROUTINE sol_mat
+
+
+   SUBROUTINE sol_exd( pt3d, cd_type )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine sol_exd  ***
+      !!                  
+      !! ** Purpose :   Reorder gcb coefficient on the extra outer  halo 
+      !!                at north fold in case of T or F pivot
+      !!
+      !! ** Method  :   Perform a circular permutation of the coefficients on 
+      !!                the total area strictly above the pivot point,
+      !!                and on the semi-row of the pivot point   
+      !!                
+      !! History :
+      !!   9.0  !  05-09  (R. Benshila)  original routine
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER(len=1) , INTENT( in ) ::   &
+         cd_type       ! define the nature of pt2d array grid-points
+         !             !  = T , U , V , F , W 
+         !             !  = S : T-point, north fold treatment
+         !             !  = G : F-point, north fold treatment
+         !             !  = I : sea-ice velocity at F-point with index shift
+      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4), INTENT( inout ) ::   &
+         pt3d          ! 2D array on which the boundary condition is applied
+
+      !! * Local variables
+      INTEGER  ::   ji, jk      ! dummy loop indices
+      INTEGER  ::   iloc                ! temporary integers
+      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4) ::   &
+         ztab          ! 2D array on which the boundary condition is applied
+      !!----------------------------------------------------------------------
+
+      ztab = pt3d
+
+      ! north fold treatment
+      ! -----------------------
+  
+      SELECT CASE ( npolj )
+         
+         CASE ( 3 , 4 )   !  T pivot
+         iloc = jpiglo/2 +1 
+            
+            SELECT CASE ( cd_type )
+  
+            CASE ( 'T', 'S', 'U', 'W' )
+               DO jk =1, 4
+                  DO ji = 1-jpr2di, nlci+jpr2di
+                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
+                  ENDDO
+               ENDDO
+
+              DO jk =1, 4
+                  DO ji = nlci+jpr2di, 1-jpr2di,  -1
+                     IF( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) &
+                       & .OR. ( mi0(iloc) == jpi+1 ) ) EXIT
+                     pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
+                  ENDDO
+               ENDDO
+
+            CASE ( 'F' ,'G' , 'I', 'V' )
+               DO jk =1, 4
+                  DO ji = 1-jpr2di, nlci+jpr2di
+                     pt3d(ji,nlcj-1:nlcj+jpr2dj,jk) = ztab(ji,nlcj-1:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
+                  ENDDO
+               ENDDO
+
+            END SELECT   ! cd_type
+  
+         CASE ( 5 , 6 )                 ! F pivot
+          iloc=jpiglo/2
+
+            SELECT CASE (cd_type )
+
+            CASE ( 'T'  ,'S', 'U', 'W')
+               DO jk =1, 4
+                  DO ji = 1-jpr2di, nlci+jpr2di
+                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
+                  ENDDO
+               ENDDO
+
+            CASE ( 'F' ,'G' , 'I', 'V' )
+               DO jk =1, 4
+                  DO ji = 1-jpr2di, nlci+jpr2di
+                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
+                  ENDDO
+               ENDDO
+               DO jk =1, 4
+                  DO ji = nlci+jpr2di, 1-jpr2di,  -1
+                    IF ( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) &
+                       & .OR. ( mi0(iloc) == jpi+1 ) ) EXIT
+                    pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
+                  ENDDO
+               ENDDO
+
+            END SELECT   ! cd_type
+
+         END SELECT   ! npolj
+  
+   END SUBROUTINE sol_exd
 
 #if defined key_feti
Index: trunk/NEMO/OPA_SRC/SOL/solver.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SOL/solver.F90	(revision 311)
+++ trunk/NEMO/OPA_SRC/SOL/solver.F90	(revision 312)
@@ -152,10 +152,25 @@
       CASE ( 1 )                ! preconditioned conjugate gradient solver
          IF(lwp) WRITE(numout,*) '          a preconditioned conjugate gradient solver is used'
+         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
+            IF(lwp) WRITE(numout,cform_err)
+            IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero'
+            nstop = nstop + 1
+         ENDIF
 
       CASE ( 2 )                ! successive-over-relaxation solver
          IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver is used'
+         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
+            IF(lwp) WRITE(numout,cform_err)
+            IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero'
+            nstop = nstop + 1
+         ENDIF
 
       CASE ( 3 )                ! FETI solver
          IF(lwp) WRITE(numout,*) '          the FETI solver is used'
+         IF( jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
+            IF(lwp) WRITE(numout,cform_err)
+            IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj should be equal to zero'
+            nstop = nstop + 1
+         ENDIF
          IF( .NOT.lk_mpp ) THEN
             IF(lwp) WRITE(numout,cform_err)
@@ -170,4 +185,21 @@
          ENDIF
          
+      CASE ( 4 )                ! successive-over-relaxation solver with extra outer halo
+         IF(lwp) WRITE(numout,*) '          a successive-over-relaxation solver with extra outer halo is used'
+         IF(lwp) WRITE(numout,*) '          with jpr2di =', jpr2di, ' and  jpr2dj =', jpr2dj
+         IF( .NOT. lk_mpp .AND. jpr2di /= 0 .AND. jpr2dj /= 0 ) THEN
+            IF(lwp) WRITE(numout,cform_err)
+            IF(lwp) WRITE(numout,*) ' jpr2di and jpr2dj are not equal to zero'
+            IF(lwp) WRITE(numout,*) ' In this case this algorithm should be used only with the key_mpp_... option'
+            nstop = nstop + 1
+         ELSE
+            IF( ( ( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) .OR. ( jpni /= 1 ) ) &
+              &  .AND. ( jpr2di /= jpr2dj ) ) THEN  
+               IF(lwp) WRITE(numout,cform_err)
+               IF(lwp) WRITE(numout,*) '          jpr2di should be equal to jpr2dj'
+               nstop = nstop + 1
+            ENDIF
+         ENDIF
+
       CASE DEFAULT
          IF(lwp) WRITE(numout,cform_err)
