source: trunk/NEMO/OPA_SRC/SOL/solsor.F90 @ 3

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

   !!----------------------------------------------------------------------
   !!   sol_sor     : 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( kt, kindic )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE sol_sor  ***
      !!                 
      !! ** Purpose :   Solve the ellipic equation for the barotropic stream 
      !!      function system (default option) or the transport divergence 
      !!      system (key_dynspg_fsc = T) using a 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
      !!   8.5  !  02-08  (G. Madec)  F90: Free form
      !!        !  03-02  (C. Deltel)  Red-Black SOR <== Not done yet!!!
      !!                                                *************
      !!----------------------------------------------------------------------
      !! * Arguments
      INTEGER, INTENT(  in   ) ::   kt       ! ocean time-step
      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
      REAL(wp) ::   zgwgt                    ! temporary scalar
      !!----------------------------------------------------------------------
      
      
      ! Iterative loop 
      ! ==============
      
      IF( kt == nit000 )  gccd(:,:) = sor * gcp(:,:,2)


      DO jn = 1, nmax

         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
         
         ! boundary conditions (at each sor iteration) only cyclic b. c. are required
#if defined key_dynspg_fsc
# if defined key_mpp
            ! Mpp: export boundary values to neighbouring processors
            CALL lbc_lnk( gcx, 'S', 1. )   ! S=T with special staff ???
# else
            CALL lbc_lnk( gcx, 'T', 1. )
# endif
#else
# if defined key_mpp
            ! Mpp: export boundary values to neighbouring processors
            CALL lbc_lnk( gcx, 'G', 1. )   ! G= F with special staff ???
# else
            CALL lbc_lnk( gcx, 'F', 1. )
# endif
#endif
         
         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
         
         ! 1. Residus
         ! ----------
 
         DO jj = 2, jpjm1
            DO ji = 2, jpim1
               gcr(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)
            END DO
         END DO

         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
         
         ! 1.1 Boundary conditions (at each sor iteration) only cyclic b. c. are required
#if defined key_dynspg_fsc
#   if defined key_mpp
            ! Mpp: export boundary values to neighbouring processors
            CALL lbc_lnk( gcr, 'S', 1. )
#   else
            ! mono- or macro-tasking: W-point, >0, 2D array, no slab
            CALL lbc_lnk( gcr, 'T', 1. )
#   endif
#else
#   if defined key_mpp
            ! Mpp: export boundary values to neighbouring processors
            CALL lbc_lnk( gcr, 'G', 1. )
#   else
            ! mono- or macro-tasking: W-point, >0, 2D array, no slab
            CALL lbc_lnk( gcr, 'F', 1. )
#   endif
#endif

         ! 1.2 Successive over relaxation
         
         DO jj = 2, jpj
            DO ji = 1, jpi
               gcr(ji,jj) = gcr(ji,jj) - sor*gcp(ji,jj,1)*gcr(ji,jj-1)
            END DO
            DO ji = 2, jpi
               gcr(ji,jj) = gcr(ji,jj) - sor*gcp(ji,jj,2)*gcr(ji-1,jj)
            END DO
         END DO
         
         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
         
         ! gcx guess
         
         DO jj = 2, jpjm1
            DO ji = 1, jpi
               gcx(ji,jj)  = (gcx(ji,jj)+sor*gcr(ji,jj))*bmask(ji,jj)
            END DO
         END DO
         
         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
         
         ! boundary conditions (at each sor iteration) only cyclic b. c. are required
#if defined key_dynspg_fsc
#   if defined key_mpp
         ! Mpp: export boundary values to neighbouring processors
         CALL lbc_lnk( gcx, 'S', 1. )
#   else
         ! mono- or macro-tasking: W-point, >0, 2D array, no slab
         CALL lbc_lnk( gcx, 'T', 1. )
#   endif
#else
#   if defined key_mpp
         ! Mpp: export boundary values to neighbouring processors
         CALL lbc_lnk( gcx, 'G', 1. )
#   else
         ! mono- or macro-tasking: W-point, >0, 2D array, no slab
         CALL lbc_lnk( gcx, 'F', 1. )
#   endif
#endif
         
         ! maximal residu (old exit test on the maximum value of residus)
         ! 
         ! imax = isamax( jpi*jpj, gcr, 1 )
         
         ! avoid an out of bound in no bounds compilation
         
         ! iimax1 = mod( imax, jpi )
         ! ijmax1 = int( float(imax) / float(jpi)) + 1
         ! resmax = abs( gcr(iimax1,ijmax1) )
         
         ! relative precision
         
         rnorme = 0.
         DO jj = 1, jpj
            DO ji = 1, jpi
               zgwgt = gcdmat(ji,jj) * gcr(ji,jj)
               rnorme= rnorme + gcr(ji,jj)*zgwgt
            END DO
         END DO
         
#if defined key_mpp
         ! mpp sum over all the global domain
         CALL  mpp_sum( rnorme )
#endif
         
         ! test of convergence
         ! old test (either res<resmax or jn=nmax)
         !           IF( res < resmax .OR. jn == nmax ) THEN
         ! relative precision
         IF( rnorme < epsr .OR. jn == nmax ) THEN
            res = SQRT( rnorme )
            niter = jn
            ncut = 999
         ENDIF
         
         !****
         !     IF(lwp)WRITE(numsol,9300) jn, res, sqrt( epsr ) / eps
9300     FORMAT('          niter :',i4,' res :',e20.10,' b :',e20.10)
         !****
         
         !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
         
         ! indicator of non-convergence or explosion
         
         IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
         IF( ncut == 999 ) GOTO 999
         
         
         ! END of iterative loop
         ! =====================
         
      END DO
      
      
999   CONTINUE
      
      
      !  2. Output in gcx
      !  -----------------
      
      ! boundary conditions (est-ce necessaire? je ne crois pas!!!!)
      
#if defined key_dynspg_fsc
# if defined key_mpp
      ! Mpp: export boundary values to neighbouring processors
      CALL lbc_lnk( gcx, 'S', 1. )
# else
      IF( nperio /= 0 ) THEN
         CALL lbc_lnk( gcx, 'T', 1. )
      ENDIF
# endif
#else
# if defined key_mpp
      ! Mpp: export boundary values to neighbouring processors
      CALL lbc_lnk( gcx, 'G', 1. )
# else
      IF( nperio /= 0 ) THEN
         CALL lbc_lnk( gcx, 'F', 1. )
      ENDIF
# endif
#endif
      
   END SUBROUTINE sol_sor

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