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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SOL/solsor.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • Property svn:keywords set to Id
File size: 7.8 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      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
61      USE wrk_nemo, ONLY: ztab => wrk_2d_1
62      !!
63      INTEGER, INTENT(inout) ::   kindic   ! solver indicator, < 0 if the convergence is not reached:
64      !                                    ! the model is stopped in step (set to zero before the call of solsor)
65      !!
66      INTEGER  ::   ji, jj, jn               ! dummy loop indices
67      INTEGER  ::   ishift, icount
68      INTEGER  ::   ijmppodd, ijmppeven, ijpr2d
69      REAL(wp) ::   ztmp, zres, zres2
70      !!----------------------------------------------------------------------
71     
72      IF( wrk_in_use(2, 1) )THEN
73         CALL ctl_stop('sol_sor: requested workspace array is unavailable')
74         RETURN
75      END IF
76
77      ijmppeven = MOD( nimpp+njmpp+jpr2di+jpr2dj   , 2 )
78      ijmppodd  = MOD( nimpp+njmpp+jpr2di+jpr2dj+1 , 2 )
79      ijpr2d    = MAX( jpr2di , jpr2dj )
80      icount = 0
81      !                                                       ! ==============
82      DO jn = 1, nn_nmax                                      ! Iterative loop
83         !                                                    ! ==============
84
85         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )   ! lateral boundary conditions
86       
87         ! Residus
88         ! -------
89
90         ! Guess black update
91         DO jj = 2-jpr2dj, nlcj-1+jpr2dj
92            ishift = MOD( jj-ijmppodd-jpr2dj, 2 )
93            DO ji = 2-jpr2di+ishift, nlci-1+jpr2di, 2
94               ztmp =                  gcb(ji  ,jj  )   &
95                  &   - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
96                  &   - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
97                  &   - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
98                  &   - gcp(ji,jj,4) * gcx(ji  ,jj+1)
99               ! Estimate of the residual
100               zres = ztmp - gcx(ji,jj)
101               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
102               ! Guess update
103               gcx(ji,jj) = rn_sor * ztmp + (1-rn_sor) * gcx(ji,jj)
104            END DO
105         END DO
106         icount = icount + 1 
107 
108         IF( MOD(icount,ijpr2d+1) == 0 )   CALL lbc_lnk_e( gcx, c_solver_pt, 1. )   ! lateral boundary conditions
109
110         ! Guess red update
111         DO jj = 2-jpr2dj, nlcj-1+jpr2dj
112            ishift = MOD( jj-ijmppeven-jpr2dj, 2 )
113            DO ji = 2-jpr2di+ishift, nlci-1+jpr2di, 2
114               ztmp =                  gcb(ji  ,jj  )   &
115                  &   - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
116                  &   - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
117                  &   - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
118                  &   - gcp(ji,jj,4) * gcx(ji  ,jj+1) 
119               ! Estimate of the residual
120               zres = ztmp - gcx(ji,jj)
121               gcr(ji,jj) = zres * gcdmat(ji,jj) * zres
122               ! Guess update
123               gcx(ji,jj) = rn_sor * ztmp + (1-rn_sor) * gcx(ji,jj)
124            END DO
125         END DO
126         icount = icount + 1
127
128         ! test of convergence
129         IF ( jn > nn_nmin .AND. MOD( jn-nn_nmin, nn_nmod ) == 0 ) THEN
130
131            SELECT CASE ( nn_sol_arp )
132            CASE ( 0 )                 ! absolute precision (maximum value of the residual)
133               zres2 = MAXVAL( gcr(2:nlci-1,2:nlcj-1) )
134               IF( lk_mpp )   CALL mpp_max( zres2 )   ! max over the global domain
135               ! test of convergence
136               IF( zres2 < rn_resmax .OR. jn == nn_nmax ) THEN
137                  res = SQRT( zres2 )
138                  niter = jn
139                  ncut = 999
140               ENDIF
141            CASE ( 1 )                 ! relative precision
142               ztab = 0.
143               ztab(:,:) = gcr(2:nlci-1,2:nlcj-1)
144               rnorme = glob_sum( ztab)    ! sum over the global domain
145               ! test of convergence
146               IF( rnorme < epsr .OR. jn == nn_nmax ) THEN
147                  res = SQRT( rnorme )
148                  niter = jn
149                  ncut = 999
150               ENDIF
151            END SELECT
152         
153         !****
154         !     IF(lwp)WRITE(numsol,9300) jn, res, sqrt( epsr ) / eps
1559300     FORMAT('          niter :',i4,' res :',e20.10,' b :',e20.10)
156         !****
157         
158         ENDIF
159         ! indicator of non-convergence or explosion
160         IF( jn == nn_nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
161         IF( ncut == 999 ) GOTO 999
162         
163         !                                                 ! =====================
164      END DO                                               ! END of iterative loop
165      !                                                    ! =====================
166     
167999   CONTINUE
168     
169      !  Output in gcx
170      !  -------------
171      CALL lbc_lnk_e( gcx, c_solver_pt, 1. )    ! boundary conditions
172      !
173      IF( wrk_not_released(2, 1) )THEN
174         CALL ctl_stop('sol_sor: failed to release workspace array')
175      END IF
176      !
177   END SUBROUTINE sol_sor
178
179   !!=====================================================================
180END MODULE solsor
Note: See TracBrowser for help on using the repository browser.