source: trunk/NEMO/OPA_SRC/SOL/solpcg.F90 @ 787

MODULE solpcg
   !!======================================================================
   !!                     ***  MODULE  solfet
   !! Ocean solver :  preconditionned conjugate gradient solver
   !!=====================================================================

   !!----------------------------------------------------------------------
   !!   sol_pcg    : preconditionned conjugate gradient solver
   !!----------------------------------------------------------------------
   !! * Modules used
   USE oce             ! ocean dynamics and tracers variables
   USE dom_oce         ! ocean space and time domain variables 
   USE sol_oce         ! ocean solver variables
   USE lib_mpp         ! distributed memory computing
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
   USE in_out_manager  ! I/O manager

   IMPLICIT NONE
   PRIVATE

   !! * Routine accessibility
   PUBLIC sol_pcg              ! ???

   !! * Substitutions
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
   !! $Header$ 
   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE sol_pcg( kindic )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE sol_pcg  ***
      !!                    
      !! ** Purpose :   Solve the ellipic equation for the barotropic stream 
      !!      function system (lk_dynspg_rl=T) or the transport divergence 
      !!      system (lk_dynspg_flt=T) using a diagonal preconditionned
      !!      conjugate gradient method.
      !!      In the former case, the barotropic stream function trend has a
      !!      zero boundary condition along all coastlines (i.e. continent
      !!      as well as islands) while in the latter the boundary condition
      !!      specification is not required.
      !!
      !! ** Method  :   Diagonal preconditionned conjugate gradient method.
      !!      the algorithm is multitasked. (case of 5 points matrix)
      !!      define              pa  = q^-1 * a
      !!                        pgcb  = q^-1 * gcb
      !!                 < . ; . >_q  = ( . )^t q ( . )
      !!      where q is the preconditioning matrix = diagonal matrix of the
      !!                                              diagonal elements of a
      !!      Initialization  :
      !!         x(o) = gcx
      !!         r(o) = d(o) = pgcb - pa.x(o)
      !!         rr(o)= < r(o) , r(o) >_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)
      !!      Convergence test :
      !!         rr(n+1) / < gcb , gcb >_q   =< epsr
      !!
      !! ** Action : - niter  : solver number of iteration done
      !!             - res    : solver residu reached
      !!             - gcx()  : solution of the elliptic system
      !!
      !! References :
      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
      !!      D Azevedo et al. 1993, Computer Science Technical Report, Tennessee U.
      !!
      !! History :
      !!        !  90-10  (G. Madec)  Original code
      !!        !  91-11  (G. Madec)
      !!        !  93-04  (M. Guyon)  loops and suppress pointers
      !!        !  95-09  (M. Imbard, J. Escobar)  mpp exchange 
      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
      !!        !  96-11  (A. Weaver)  correction to preconditioning
      !!   8.5  !  02-08  (G. Madec)  F90: Free form
      !!        !  08-01  (R. Benshila) mpp optimization
      !!----------------------------------------------------------------------
      !! * Arguments
      INTEGER, INTENT( inout ) ::   kindic   ! solver indicator, < 0 if the conver-
      !                                      ! gence is not reached: the model is
      !                                      ! stopped in step
      !                                      ! set to zero before the call of solpcg

      !! * Local declarations
      INTEGER ::   ji, jj, jn                ! dummy loop indices
      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
         !                                             !================

         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.
               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
 
         ! rnorme = (gcr,gcr)
         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
            res = SQRT( rnorme )
            niter = jn
            ncut = 999
         ENDIF

         ! beta = (rk+1,rk+1)/(rk,rk)
         beta = rnorme / rr
         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

         !                                             !================
      END DO                                           !    End Loop
      !                                                !================
     
999   CONTINUE
     
     
      ! Output in gcx with lateral b.c. applied
      ! ---------------------------------------
     
      CALL lbc_lnk( gcx, c_solver_pt, 1. )
     
   END SUBROUTINE sol_pcg

   !!=====================================================================
END MODULE solpcg