Index: /trunk/NEMO/OPA_SRC/SOL/solpcg.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/SOL/solpcg.F90	(revision 786)
+++ /trunk/NEMO/OPA_SRC/SOL/solpcg.F90	(revision 787)
@@ -52,16 +52,21 @@
       !!      where q is the preconditioning matrix = diagonal matrix of the
       !!                                              diagonal elements of a
-      !!      Initialization:
+      !!      Initialization  :
       !!         x(o) = gcx
       !!         r(o) = d(o) = pgcb - pa.x(o)
       !!         rr(o)= < r(o) , r(o) >_q
-      !!      Iteration n   :
-      !!         z(n)   = pa.d(n)
-      !!         alp(n) = rr(n) / < z(n) , d(n) >_q
+      !!      Iteration 1     :
+      !!         standard PCG algorithm
+      !!      Iteration n > 1 :
+      !!         s(n)   = pa.r(n)
+      !!         gam(n) = < r(n) , r(n) >_q
+      !!         del(n) = < r(n) , s(n) >_q
+      !!         bet(n) = gam(n) / gam(n-1)
+      !!         d(n)   = r(n) + bet(n) d(n-1)
+      !!         z(n)   = s(n) + bet(n) z(n-1) 
+      !!         sig(n) = del(n) - bet(n)*bet(n)*sig(n-1) 
+      !!         alp(n) = gam(n) / sig(n) 
       !!         x(n+1) = x(n) + alp(n) d(n)
       !!         r(n+1) = r(n) - alp(n) z(n)
-      !!         rr(n+1)= < r(n+1) , r(n+1) >_q
-      !!         bet(n) = rr(n+1) / rr(n)
-      !!         r(n+1) = r(n+1) + bet(n+1) d(n)
       !!      Convergence test :
       !!         rr(n+1) / < gcb , gcb >_q   =< epsr
@@ -73,4 +78,5 @@
       !! References :
       !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
+      !!      D Azevedo et al. 1993, Computer Science Technical Report, Tennessee U.
       !!
       !! History :
@@ -82,4 +88,5 @@
       !!        !  96-11  (A. Weaver)  correction to preconditioning
       !!   8.5  !  02-08  (G. Madec)  F90: Free form
+      !!        !  08-01  (R. Benshila) mpp optimization
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -91,66 +98,87 @@
       !! * Local declarations
       INTEGER ::   ji, jj, jn                ! dummy loop indices
-      REAL(wp) ::   zgcad                    ! temporary scalars
+      REAL(wp) ::  zgcad                     ! temporary scalars
+      REAL(wp), DIMENSION(2) :: zsum
+      REAL(wp), DIMENSION(jpi,jpj) :: zgcr
       !!----------------------------------------------------------------------
 
+      ! Initialization of the algorithm with standard PCG
+      ! -------------------------------------------------
+
+      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
+
+      ! gcr   = gcb-a.gcx
+      ! gcdes = gcr
+
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
+               &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
+               &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
+               &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
+               &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   )
+            gcr  (ji,jj) = zgcad
+            gcdes(ji,jj) = zgcad
+         END DO
+      END DO
+
+      ! rnorme = (gcr,gcr)
+      rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
+      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
+
+      CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
+
+      ! gccd = matrix . gcdes
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   &
+               &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
+               &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   )
+         END DO
+      END DO 
+
+      ! alph = (gcr,gcr)/(gcdes,gccd)
+      radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
+      IF( lk_mpp )   CALL mpp_sum( radd )   ! sum over the global domain
+      alph = rnorme /radd
+
+      ! gcx = gcx + alph * gcdes
+      ! gcr = gcr - alph * gccd
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
+            gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
+         END DO
+      END DO
+
+      ! Algorithm wtih Eijkhout rearrangement
+      ! -------------------------------------
+        
       !                                                !================
       DO jn = 1, nmax                                  ! Iterative loop
          !                                             !================
 
-         IF( jn == 1 ) THEN           ! Initialization of the algorithm
-
-            CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
-
-            ! gcr   = gcb-a.gcx
-            ! gcdes = gsr
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
-                     &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
-                     &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
-                     &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
-                     &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   )
-                  gcr  (ji,jj) = zgcad
-                  gcdes(ji,jj) = zgcad
-               END DO
-            END DO
-            
-            rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
-            IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
-            rr = rnorme
-
-         ENDIF
-        
-         !                             ! Algorithm
-        
-         CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
-        
-         ! ... gccd = matrix . gcdes
+         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! lateral boundary condition
+
+         ! zgcr = matrix . gcr
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   &
-                  &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
-                  &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   )
+               zgcr(ji,jj) = bmask(ji,jj)*( gcr(ji,jj)   &
+                  &        +gcp(ji,jj,1)*gcr(ji,jj-1)+gcp(ji,jj,2)*gcr(ji-1,jj)   &
+                  &        +gcp(ji,jj,4)*gcr(ji,jj+1)+gcp(ji,jj,3)*gcr(ji+1,jj)   )
             END DO
          END DO
  
-         ! alph = (gcr,gcr)/(gcdes,gccd)
-         radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
-         IF( lk_mpp )   CALL mpp_sum( radd )   ! sum over the global domain
-         alph = rr / radd
-         
-         ! gcx = gcx + alph * gcdes
-         ! gcr = gcr - alph * gccd
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
-               gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
-            END DO
-         END DO
-        
          ! rnorme = (gcr,gcr)
-         rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
-         IF( lk_mpp )   CALL  mpp_sum( rnorme )   ! sum over the global domain
-        
+         rr = rnorme
+         zsum(1) = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
+
+         ! zgcad = (zgcr,gcr) 
+         zsum(2) = SUM( gcr(2:jpim1,2:jpjm1) * gcdmat(2:jpim1,2:jpjm1) * zgcr(2:jpim1,2:jpjm1) )
+
+         IF( lk_mpp )   CALL mpp_sum( zsum, 2 )   ! sum over the global domain
+         rnorme = zsum(1)  
+         zgcad  = zsum(2)
+
          ! test of convergence
          IF( rnorme < epsr .OR. jn == nmax ) THEN
@@ -159,20 +187,25 @@
             ncut = 999
          ENDIF
-        
+
          ! beta = (rk+1,rk+1)/(rk,rk)
          beta = rnorme / rr
-         rr   = rnorme
-
+         radd = zgcad - beta*beta*radd
+         alph = rnorme / radd
+
+         ! gcx = gcx + alph * gcdes
+         ! gcr = gcr - alph * gccd
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               gcdes(ji,jj) = gcr (ji,jj) + beta * gcdes(ji,jj) 
+               gccd (ji,jj) = zgcr(ji,jj) + beta * gccd (ji,jj) 
+               gcx  (ji,jj) = gcx (ji,jj) + alph * gcdes(ji,jj) 
+               gcr  (ji,jj) = gcr (ji,jj) - alph * gccd (ji,jj) 
+            END DO
+         END DO
+        
          ! indicator of non-convergence or explosion
          IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
          IF( ncut == 999 ) GOTO 999
 
-         ! gcdes = gcr + beta * gcdes
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               gcdes(ji,jj) = bmask(ji,jj)*( gcr(ji,jj) + beta * gcdes(ji,jj) )
-            END DO
-         END DO
-        
          !                                             !================
       END DO                                           !    End Loop
