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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SOL/solsor.F90 @ 4460

Last change on this file since 4460 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 8.0 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   !! * Control permutation of array indices
32#  include "oce_ftrans.h90"
33#  include "dom_oce_ftrans.h90"
34#  include "zdf_oce_ftrans.h90"
35
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42     
43   SUBROUTINE sol_sor( kindic )
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE sol_sor  ***
46      !!                 
47      !! ** Purpose :   Solve the ellipic equation for the transport
48      !!      divergence system  using a red-black successive-over-
49      !!      relaxation method.
50      !!       This routine provides a MPI optimization to the existing solsor
51      !!     by reducing the number of call to lbc.
52      !!
53      !! ** Method  :   Successive-over-relaxation method using the red-black
54      !!      technique. The former technique used was not compatible with
55      !!      the north-fold boundary condition used in orca configurations.
56      !!      Compared to the classical sol_sor, this routine provides a
57      !!      mpp optimization by reducing the number of calls to lnc_lnk
58      !!      The solution is computed on a larger area and the boudary
59      !!      conditions only when the inside domain is reached.
60      !!
61      !! References :   Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
62      !!                Beare and Stevens 1997 Ann. Geophysicae 15, 1369-1377
63      !!----------------------------------------------------------------------
64      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
65      USE wrk_nemo, ONLY:   ztab => wrk_2d_1    ! 2D workspace
66      !!
67      INTEGER, INTENT(inout) ::   kindic   ! solver indicator, < 0 if the convergence is not reached:
68      !                                    ! the model is stopped in step (set to zero before the call of solsor)
69      !!
70      INTEGER  ::   ji, jj, jn       ! dummy loop indices
71      INTEGER  ::   ishift, icount, ijmppodd, ijmppeven, ijpr2d   ! local integers
72      REAL(wp) ::   ztmp, zres, zres2                             ! local scalars
73      !!----------------------------------------------------------------------
74     
75      IF( wrk_in_use(2, 1) )THEN
76         CALL ctl_stop('sol_sor: requested workspace array is unavailable')   ;   RETURN
77      ENDIF
78
79      ijmppeven = MOD( nimpp+njmpp+jpr2di+jpr2dj   , 2 )
80      ijmppodd  = MOD( nimpp+njmpp+jpr2di+jpr2dj+1 , 2 )
81      ijpr2d    = MAX( jpr2di , jpr2dj )
82      icount = 0
83      !                                                       ! ==============
84      DO jn = 1, nn_nmax                                      ! Iterative loop
85         !                                                    ! ==============
86
87         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )   ! lateral boundary conditions
88       
89         ! Residus
90         ! -------
91
92         ! Guess black update
93         DO jj = 2-jpr2dj, nlcj-1+jpr2dj
94            ishift = MOD( jj-ijmppodd-jpr2dj, 2 )
95            DO ji = 2-jpr2di+ishift, nlci-1+jpr2di, 2
96               ztmp =                  gcb(ji  ,jj  )   &
97                  &   - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
98                  &   - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
99                  &   - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
100                  &   - gcp(ji,jj,4) * gcx(ji  ,jj+1)
101               ! Estimate of the residual
102               zres = ztmp - gcx(ji,jj)
103               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
104               ! Guess update
105               gcx(ji,jj) = rn_sor * ztmp + (1-rn_sor) * gcx(ji,jj)
106            END DO
107         END DO
108         icount = icount + 1 
109 
110         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )   ! lateral boundary conditions
111
112         ! Guess red update
113         DO jj = 2-jpr2dj, nlcj-1+jpr2dj
114            ishift = MOD( jj-ijmppeven-jpr2dj, 2 )
115            DO ji = 2-jpr2di+ishift, nlci-1+jpr2di, 2
116               ztmp =                  gcb(ji  ,jj  )   &
117                  &   - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
118                  &   - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
119                  &   - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
120                  &   - gcp(ji,jj,4) * gcx(ji  ,jj+1) 
121               ! Estimate of the residual
122               zres = ztmp - gcx(ji,jj)
123               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
124               ! Guess update
125               gcx(ji,jj) = rn_sor * ztmp + (1-rn_sor) * gcx(ji,jj)
126            END DO
127         END DO
128         icount = icount + 1
129
130         ! test of convergence
131         IF ( jn > nn_nmin .AND. MOD( jn-nn_nmin, nn_nmod ) == 0 ) THEN
132
133            SELECT CASE ( nn_sol_arp )
134            CASE ( 0 )                 ! absolute precision (maximum value of the residual)
135               zres2 = MAXVAL( gcr(2:nlci-1,2:nlcj-1) )
136               IF( lk_mpp )   CALL mpp_max( zres2 )   ! max over the global domain
137               ! test of convergence
138               IF( zres2 < rn_resmax .OR. jn == nn_nmax ) THEN
139                  res = SQRT( zres2 )
140                  niter = jn
141                  ncut = 999
142               ENDIF
143            CASE ( 1 )                 ! relative precision
144               ztab = 0.
145               ztab(:,:) = gcr(2:nlci-1,2:nlcj-1)
146               rnorme = glob_sum( ztab)    ! sum over the global domain
147               ! test of convergence
148               IF( rnorme < epsr .OR. jn == nn_nmax ) THEN
149                  res = SQRT( rnorme )
150                  niter = jn
151                  ncut = 999
152               ENDIF
153            END SELECT
154         
155         !****
156         !     IF(lwp)WRITE(numsol,9300) jn, res, sqrt( epsr ) / eps
1579300     FORMAT('          niter :',i4,' res :',e20.10,' b :',e20.10)
158         !****
159         
160         ENDIF
161         ! indicator of non-convergence or explosion
162         IF( jn == nn_nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
163         IF( ncut == 999 ) GOTO 999
164         
165         !                                                 ! =====================
166      END DO                                               ! END of iterative loop
167      !                                                    ! =====================
168     
169999   CONTINUE
170     
171      !  Output in gcx
172      !  -------------
173      CALL lbc_lnk_e( gcx, c_solver_pt, 1. )    ! boundary conditions
174      !
175      IF( wrk_not_released(2, 1) )   CALL ctl_stop('sol_sor: failed to release workspace array')
176      !
177   END SUBROUTINE sol_sor
178
179   !!=====================================================================
180END MODULE solsor
Note: See TracBrowser for help on using the repository browser.