source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/SOL/solsor.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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