source: tags/nemo_dev_x2/NEMO/OPA_SRC/SOL/solsor.F90 @ 9261

MODULE solsor
   !!======================================================================
   !!                     ***  MODULE  solsor  ***
   !! Ocean solver :  Successive Over-Relaxation solver
   !!=====================================================================

   !!----------------------------------------------------------------------
   !!   sol_sor     : Red-Black Successive Over-Relaxation solver
   !!----------------------------------------------------------------------
   !! * Modules used
   USE oce             ! ocean dynamics and tracers variables
   USE dom_oce         ! ocean space and time domain variables 
   USE zdf_oce         ! ocean vertical physics variables
   USE sol_oce         ! solver variables
   USE in_out_manager  ! I/O manager
   USE lib_mpp         ! distributed memory computing
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)

   IMPLICIT NONE
   PRIVATE

   !! * Routine accessibility
   PUBLIC sol_sor              ! ???

   !!----------------------------------------------------------------------
   !!   OPA 9.0 , LODYC-IPSL  (2003)
   !!----------------------------------------------------------------------

CONTAINS
      
   SUBROUTINE sol_sor( kindic )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE sol_sor  ***
      !!                 
      !! ** Purpose :   Solve the ellipic equation for the barotropic stream 
      !!      function system (lk_dynspg_rl=T) or the transport divergence 
      !!      system (lk_dynspg_fsc=T) using a red-black successive-over-
      !!      relaxation 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  :   Successive-over-relaxation method using the red-black 
      !!      technique. The former technique used was not compatible with 
      !!      the north-fold boundary condition used in orca configurations.
      !!
      !! References :
      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
      !!
      !! History :
      !!        !  90-10  (G. Madec)  Original code
      !!        !  91-11  (G. Madec)
      !!   7.1  !  93-04  (G. Madec)  time filter
      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
      !!        !  96-11  (A. Weaver)  correction to preconditioning
      !!   9.0  !  03-04  (C. Deltel, G. Madec)  Red-Black SOR in free form
      !!----------------------------------------------------------------------
      !! * 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 solsor
      !! * Local declarations
      INTEGER  ::   ji, jj, jn               ! dummy loop indices
      INTEGER  ::   ishift
      REAL(wp) ::   ztmp, zres, zres2

      INTEGER  ::   ijmppodd, ijmppeven
      !!----------------------------------------------------------------------
      
      ijmppeven = MOD(njmpp  ,2)
      ijmppodd  = MOD(njmpp+1,2)
      !                                                       ! ==============
      DO jn = 1, nmax                                         ! Iterative loop 
         !                                                    ! ==============

         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions
         
         ! Residus
         ! -------

         zres2 = 0.e0

         ! Guess black update
        DO jj = 2, jpjm1
            ishift = MOD( jj-ijmppodd, 2 )
            DO ji = 2+ishift, jpim1, 2
               ztmp =                  gcb(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)
               ! Estimate of the residual
               zres  = ztmp - gcx(ji,jj)
               zres2 = zres2 + zres * gcdmat(ji,jj) * zres
               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
               ! Guess update
               gcx(ji,jj) = sor * ztmp + (1-sor) * gcx(ji,jj)
            END DO
         END DO

         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions

         ! Guess red update
         DO jj = 2, jpjm1
            ishift = MOD( jj-ijmppeven, 2 )
            DO ji = 2+ishift, jpim1, 2
               ztmp =                  gcb(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) 
               ! Estimate of the residual
               zres  = ztmp - gcx(ji,jj)
               zres2 = zres2 + zres * gcdmat(ji,jj) * zres
               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
               ! Guess update
               gcx(ji,jj) = sor * ztmp  +  (1-sor) * gcx(ji,jj)
            END DO
         END DO

         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions

         ! relative precision
         rnorme = zres2
!!!!     IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
!i
         zres2 = MAXVAL( gcr(2:jpim1,2:jpjm1) )

         IF( lk_mpp )   CALL mpp_max( zres2 )   ! sum over the global domain
!i
         ! test of convergence
         ! old test (either res<resmax or jn=nmax)
         !           IF( res < resmax .OR. jn == nmax ) THEN
                     IF( zres2 < 1.e-10 .OR. jn == nmax ) THEN
         ! relative precision
         !        IF( rnorme < epsr .OR. jn == nmax ) THEN
         !           res = SQRT( rnorme )
            res = SQRT( zres2 )
            niter = jn
            ncut = 999
            write(numout,*) 'solsor : res max = ', zres2, 'niter= ', niter
         ENDIF
         
         !****
         !     IF(lwp)WRITE(numsol,9300) jn, res, sqrt( epsr ) / eps
9300     FORMAT('          niter :',i4,' res :',e20.10,' b :',e20.10)
         !****
         
         ! 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 of iterative loop
      !                                                    ! =====================
      
999   CONTINUE
      
      
      !  Output in gcx
      !  -------------
      
      CALL lbc_lnk( gcx, c_solver_pt, 1. )    ! boundary conditions (est-ce necessaire? je ne crois pas!!!!)
      
   END SUBROUTINE sol_sor

   !!=====================================================================
END MODULE solsor