Changeset 111 for trunk/NEMO/OPA_SRC/SOL/solsor.F90

Timestamp:

2004-06-28T15:22:55+02:00 (20 years ago)

Author:

opalod

Message:

CT : UPDATE070 : Optimisation of the Red-Black SOR algorithm convergence

File:

• trunk/NEMO/OPA_SRC/SOL/solsor.F90 (modified) (5 diffs)

trunk/NEMO/OPA_SRC/SOL/solsor.F90

@@ -76 +76 @@
          !                                                    ! ==============
 
-         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions
+         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boundary conditions
          
          ! Residus
          ! -------
 
-         zres2 = 0.e0
-
          ! Guess black update
-        DO jj = 2, jpjm1
+         DO jj = 2, jpjm1
             ishift = MOD( jj-ijmppodd, 2 )
             DO ji = 2+ishift, jpim1, 2
@@ -93 +91 @@
                   &   - gcp(ji,jj,4) * gcx(ji  ,jj+1)
                ! Estimate of the residual
-               zres  = ztmp - gcx(ji,jj)
-               zres2 = zres2 + zres * gcdmat(ji,jj) * zres
+               zres = ztmp - gcx(ji,jj)
                gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
                ! Guess update
@@ -113 +110 @@
                   &   - gcp(ji,jj,4) * gcx(ji  ,jj+1) 
                ! Estimate of the residual
-               zres  = ztmp - gcx(ji,jj)
-               zres2 = zres2 + zres * gcdmat(ji,jj) * zres
+               zres = ztmp - gcx(ji,jj)
                gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
                ! Guess update
-               gcx(ji,jj) = sor * ztmp  +  (1-sor) * gcx(ji,jj)
+               gcx(ji,jj) = sor * ztmp + (1-sor) * gcx(ji,jj)
             END DO
          END DO
@@ -123 +119 @@
          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) )
+         ! test of convergence
+         IF ( jn > nmin .AND. MOD( jn-nmin, nmod ) == 0 ) then
 
-         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
+            SELECT CASE ( nsol_arp )
+            CASE ( 1 )                 ! relative precision
+               rnorme = SUM( gcr(2:jpim1,2:jpjm1) )
+               IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
+               ! test of convergence
+               IF( rnorme < epsr .OR. jn == nmax ) THEN
+                  res = SQRT( rnorme )
+                  niter = jn
+                  ncut = 999
+               ENDIF
+            CASE ( 0 )                 ! absolute precision (maximum value of the residual)
+               zres2 = MAXVAL( gcr(2:jpim1,2:jpjm1) )
+               IF( lk_mpp )   CALL mpp_max( zres2 )   ! max over the global domain
+               ! test of convergence
+               IF( zres2 < resmax .OR. jn == nmax ) THEN
+                  res = SQRT( zres2 )
+                  niter = jn
+                  ncut = 999
+               ENDIF
+            END CASE
          
          !****
@@ -149 +148 @@
          !****
          
+         ENDIF
          ! indicator of non-convergence or explosion
          IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2