Changeset 86 for trunk/NEMO/OPA_SRC/SOL

Timestamp:

2004-04-22T15:37:51+02:00 (20 years ago)

Author:

opalod

Message:

CT : BUGFIX058 : Red-Black Successive Over-Relaxation solver runs in mpp

File:

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

trunk/NEMO/OPA_SRC/SOL/solsor.F90

@@ -1 +1 @@
 MODULE solsor
    !!======================================================================
-   !!                     ***  MODULE  solsor
+   !!                     ***  MODULE  solsor  ***
    !! Ocean solver :  Successive Over-Relaxation solver
    !!=====================================================================
 
    !!----------------------------------------------------------------------
-   !!   sol_sor     : Successive Over-Relaxation solver
+   !!   sol_sor     : Red-Black Successive Over-Relaxation solver
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -35 +35 @@
       !! ** 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 successive-over-relaxation
-      !!      method.
+      !!      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
@@ -55 +55 @@
       !!        !  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!!!
-      !!                                                *************
+      !!   9.0  !  03-04  (C. Deltel, G. Madec)  Red-Black SOR in free form
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -64 +62 @@
       !                                      ! 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 
@@ -77 +80 @@
          ! 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)
+
+         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( gcr, c_solver_pt, 1. )   ! applied the lateral boubary conditions
 
-         ! 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)
+         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
-         
-         ! 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
-         CALL lbc_lnk( gcx, c_solver_pt, 1. )
-         
+
+         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions
+
          ! relative precision
-         rnorme = SUM( gcr(:,:) * gcdmat(:,:) * gcr(:,:) )
-         IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
-         
+         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 )
+         !        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