New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
solsor.F90 in trunk/NEMO/OPA_SRC/SOL – NEMO

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

Last change on this file since 95 was 95, checked in by opalod, 19 years ago

CT : BUGFIX061 : mpp correction on ijmppeven and ijmppodd indices

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.0 KB
Line 
1MODULE solsor
2   !!======================================================================
3   !!                     ***  MODULE  solsor  ***
4   !! Ocean solver :  Successive Over-Relaxation solver
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   sol_sor     : Red-Black Successive Over-Relaxation solver
9   !!----------------------------------------------------------------------
10   !! * Modules used
11   USE oce             ! ocean dynamics and tracers variables
12   USE dom_oce         ! ocean space and time domain variables
13   USE zdf_oce         ! ocean vertical physics variables
14   USE sol_oce         ! solver variables
15   USE in_out_manager  ! I/O manager
16   USE lib_mpp         ! distributed memory computing
17   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Routine accessibility
23   PUBLIC sol_sor              ! ???
24
25   !!----------------------------------------------------------------------
26   !!   OPA 9.0 , LODYC-IPSL  (2003)
27   !!----------------------------------------------------------------------
28
29CONTAINS
30     
31   SUBROUTINE sol_sor( kindic )
32      !!----------------------------------------------------------------------
33      !!                  ***  ROUTINE sol_sor  ***
34      !!                 
35      !! ** Purpose :   Solve the ellipic equation for the barotropic stream
36      !!      function system (lk_dynspg_rl=T) or the transport divergence
37      !!      system (lk_dynspg_fsc=T) using a red-black successive-over-
38      !!      relaxation method.
39      !!       In the former case, the barotropic stream function trend has a
40      !!     zero boundary condition along all coastlines (i.e. continent
41      !!     as well as islands) while in the latter the boundary condition
42      !!     specification is not required.
43      !!
44      !! ** Method  :   Successive-over-relaxation method using the red-black
45      !!      technique. The former technique used was not compatible with
46      !!      the north-fold boundary condition used in orca configurations.
47      !!
48      !! References :
49      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
50      !!
51      !! History :
52      !!        !  90-10  (G. Madec)  Original code
53      !!        !  91-11  (G. Madec)
54      !!   7.1  !  93-04  (G. Madec)  time filter
55      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
56      !!        !  96-11  (A. Weaver)  correction to preconditioning
57      !!   9.0  !  03-04  (C. Deltel, G. Madec)  Red-Black SOR in free form
58      !!----------------------------------------------------------------------
59      !! * Arguments
60      INTEGER, INTENT( inout ) ::   kindic   ! solver indicator, < 0 if the conver-
61      !                                      ! gence is not reached: the model is
62      !                                      ! stopped in step
63      !                                      ! set to zero before the call of solsor
64      !! * Local declarations
65      INTEGER  ::   ji, jj, jn               ! dummy loop indices
66      INTEGER  ::   ishift
67      REAL(wp) ::   ztmp, zres, zres2
68
69      INTEGER  ::   ijmppodd, ijmppeven
70      !!----------------------------------------------------------------------
71     
72      ijmppeven = MOD(nimpp+njmpp  ,2)
73      ijmppodd  = MOD(nimpp+njmpp+1,2)
74      !                                                       ! ==============
75      DO jn = 1, nmax                                         ! Iterative loop
76         !                                                    ! ==============
77
78         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions
79         
80         ! Residus
81         ! -------
82
83         zres2 = 0.e0
84
85         ! Guess black update
86        DO jj = 2, jpjm1
87            ishift = MOD( jj-ijmppodd, 2 )
88            DO ji = 2+ishift, jpim1, 2
89               ztmp =                  gcb(ji  ,jj  )   &
90                  &   - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
91                  &   - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
92                  &   - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
93                  &   - gcp(ji,jj,4) * gcx(ji  ,jj+1)
94               ! Estimate of the residual
95               zres  = ztmp - gcx(ji,jj)
96               zres2 = zres2 + zres * gcdmat(ji,jj) * zres
97               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
98               ! Guess update
99               gcx(ji,jj) = sor * ztmp + (1-sor) * gcx(ji,jj)
100            END DO
101         END DO
102
103         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions
104
105         ! Guess red update
106         DO jj = 2, jpjm1
107            ishift = MOD( jj-ijmppeven, 2 )
108            DO ji = 2+ishift, jpim1, 2
109               ztmp =                  gcb(ji  ,jj  )   &
110                  &   - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
111                  &   - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
112                  &   - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
113                  &   - gcp(ji,jj,4) * gcx(ji  ,jj+1) 
114               ! Estimate of the residual
115               zres  = ztmp - gcx(ji,jj)
116               zres2 = zres2 + zres * gcdmat(ji,jj) * zres
117               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
118               ! Guess update
119               gcx(ji,jj) = sor * ztmp  +  (1-sor) * gcx(ji,jj)
120            END DO
121         END DO
122
123         CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! applied the lateral boubary conditions
124
125         ! relative precision
126         rnorme = zres2
127!!!!     IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
128!i
129         zres2 = MAXVAL( gcr(2:jpim1,2:jpjm1) )
130
131         IF( lk_mpp )   CALL mpp_max( zres2 )   ! sum over the global domain
132!i
133         ! test of convergence
134         ! old test (either res<resmax or jn=nmax)
135         !           IF( res < resmax .OR. jn == nmax ) THEN
136                     IF( zres2 < 1.e-10 .OR. jn == nmax ) THEN
137         ! relative precision
138         !        IF( rnorme < epsr .OR. jn == nmax ) THEN
139         !           res = SQRT( rnorme )
140            res = SQRT( zres2 )
141            niter = jn
142            ncut = 999
143         !  write(numout,*) 'solsor : res max = ', zres2, 'niter= ', niter
144         ENDIF
145         
146         !****
147         !     IF(lwp)WRITE(numsol,9300) jn, res, sqrt( epsr ) / eps
1489300     FORMAT('          niter :',i4,' res :',e20.10,' b :',e20.10)
149         !****
150         
151         ! indicator of non-convergence or explosion
152         IF( jn == nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
153         IF( ncut == 999 ) GOTO 999
154         
155         !                                                 ! =====================
156      END DO                                               ! END of iterative loop
157      !                                                    ! =====================
158     
159999   CONTINUE
160     
161     
162      !  Output in gcx
163      !  -------------
164     
165      CALL lbc_lnk( gcx, c_solver_pt, 1. )    ! boundary conditions (est-ce necessaire? je ne crois pas!!!!)
166     
167   END SUBROUTINE sol_sor
168
169   !!=====================================================================
170END MODULE solsor
Note: See TracBrowser for help on using the repository browser.