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/NEMOGCM/NEMO/OPA_SRC/SOL – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/SOL/solsor.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

  • Property svn:keywords set to Id
File size: 7.5 KB
Line 
1MODULE solsor
2   !!======================================================================
3   !!                     ***  MODULE  solsor  ***
4   !! Ocean solver :  Successive Over-Relaxation solver
5   !!=====================================================================
6   !! History :  OPA  ! 1990-10  (G. Madec)  Original code
7   !!            7.1  ! 1993-04  (G. Madec)  time filter
8   !!                 ! 1996-05  (G. Madec)  merge sor and pcg formulations
9   !!                 ! 1996-11  (A. Weaver)  correction to preconditioning
10   !!   NEMO     1.0  ! 2003-04  (C. Deltel, G. Madec)  Red-Black SOR in free form
11   !!            2.0  ! 2005-09  (R. Benshila, G. Madec)  MPI optimization
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   sol_sor     : Red-Black Successive Over-Relaxation solver
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers variables
18   USE dom_oce         ! ocean space and time domain variables
19   USE zdf_oce         ! ocean vertical physics variables
20   USE sol_oce         ! solver variables
21   USE in_out_manager  ! I/O manager
22   USE lib_mpp         ! distributed memory computing
23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
24   USE lib_fortran
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   sol_sor    !
30
31   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
33   !! $Id$
34   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36
37CONTAINS
38     
39   SUBROUTINE sol_sor( kindic )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE sol_sor  ***
42      !!                 
43      !! ** Purpose :   Solve the ellipic equation for the transport
44      !!      divergence system  using a red-black successive-over-
45      !!      relaxation method.
46      !!       This routine provides a MPI optimization to the existing solsor
47      !!     by reducing the number of call to lbc.
48      !!
49      !! ** Method  :   Successive-over-relaxation method using the red-black
50      !!      technique. The former technique used was not compatible with
51      !!      the north-fold boundary condition used in orca configurations.
52      !!      Compared to the classical sol_sor, this routine provides a
53      !!      mpp optimization by reducing the number of calls to lnc_lnk
54      !!      The solution is computed on a larger area and the boudary
55      !!      conditions only when the inside domain is reached.
56      !!
57      !! References :   Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
58      !!                Beare and Stevens 1997 Ann. Geophysicae 15, 1369-1377
59      !!----------------------------------------------------------------------
60      INTEGER, INTENT(inout) ::   kindic   ! solver indicator, < 0 if the convergence is not reached:
61      !                                    ! the model is stopped in step (set to zero before the call of solsor)
62      !!
63      INTEGER  ::   ji, jj, jn               ! dummy loop indices
64      INTEGER  ::   ishift, icount
65      INTEGER  ::   ijmppodd, ijmppeven, ijpr2d
66      REAL(wp) ::   ztmp, zres, zres2
67      REAL(wp), DIMENSION(jpi,jpj) ::ztab
68      !!----------------------------------------------------------------------
69     
70      ijmppeven = MOD( nimpp+njmpp+jpr2di+jpr2dj   , 2 )
71      ijmppodd  = MOD( nimpp+njmpp+jpr2di+jpr2dj+1 , 2 )
72      ijpr2d    = MAX( jpr2di , jpr2dj )
73      icount = 0
74      !                                                       ! ==============
75      DO jn = 1, nn_nmax                                      ! Iterative loop
76         !                                                    ! ==============
77
78         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )   ! lateral boundary conditions
79       
80         ! Residus
81         ! -------
82
83         ! Guess black update
84         DO jj = 2-jpr2dj, nlcj-1+jpr2dj
85            ishift = MOD( jj-ijmppodd-jpr2dj, 2 )
86            DO ji = 2-jpr2di+ishift, nlci-1+jpr2di, 2
87               ztmp =                  gcb(ji  ,jj  )   &
88                  &   - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
89                  &   - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
90                  &   - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
91                  &   - gcp(ji,jj,4) * gcx(ji  ,jj+1)
92               ! Estimate of the residual
93               zres = ztmp - gcx(ji,jj)
94               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
95               ! Guess update
96               gcx(ji,jj) = rn_sor * ztmp + (1-rn_sor) * gcx(ji,jj)
97            END DO
98         END DO
99         icount = icount + 1 
100 
101         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )   ! lateral boundary conditions
102
103         ! Guess red update
104         DO jj = 2-jpr2dj, nlcj-1+jpr2dj
105            ishift = MOD( jj-ijmppeven-jpr2dj, 2 )
106            DO ji = 2-jpr2di+ishift, nlci-1+jpr2di, 2
107               ztmp =                  gcb(ji  ,jj  )   &
108                  &   - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
109                  &   - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
110                  &   - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
111                  &   - gcp(ji,jj,4) * gcx(ji  ,jj+1) 
112               ! Estimate of the residual
113               zres = ztmp - gcx(ji,jj)
114               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
115               ! Guess update
116               gcx(ji,jj) = rn_sor * ztmp + (1-rn_sor) * gcx(ji,jj)
117            END DO
118         END DO
119         icount = icount + 1
120
121         ! test of convergence
122         IF ( jn > nn_nmin .AND. MOD( jn-nn_nmin, nn_nmod ) == 0 ) THEN
123
124            SELECT CASE ( nn_sol_arp )
125            CASE ( 0 )                 ! absolute precision (maximum value of the residual)
126               zres2 = MAXVAL( gcr(2:nlci-1,2:nlcj-1) )
127               IF( lk_mpp )   CALL mpp_max( zres2 )   ! max over the global domain
128               ! test of convergence
129               IF( zres2 < rn_resmax .OR. jn == nn_nmax ) THEN
130                  res = SQRT( zres2 )
131                  niter = jn
132                  ncut = 999
133               ENDIF
134            CASE ( 1 )                 ! relative precision
135               ztab = 0.
136               ztab(:,:) = gcr(2:nlci-1,2:nlcj-1)
137               rnorme = glob_sum( ztab)    ! sum over the global domain
138               ! test of convergence
139               IF( rnorme < epsr .OR. jn == nn_nmax ) THEN
140                  res = SQRT( rnorme )
141                  niter = jn
142                  ncut = 999
143               ENDIF
144            END SELECT
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         ENDIF
152         ! indicator of non-convergence or explosion
153         IF( jn == nn_nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
154         IF( ncut == 999 ) GOTO 999
155         
156         !                                                 ! =====================
157      END DO                                               ! END of iterative loop
158      !                                                    ! =====================
159     
160999   CONTINUE
161     
162      !  Output in gcx
163      !  -------------
164      CALL lbc_lnk_e( gcx, c_solver_pt, 1. )    ! boundary conditions
165      !
166   END SUBROUTINE sol_sor
167
168   !!=====================================================================
169END MODULE solsor
Note: See TracBrowser for help on using the repository browser.