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.
solpcg.F90 in trunk/NEMO/OPA_SRC/SOL – NEMO

source: trunk/NEMO/OPA_SRC/SOL/solpcg.F90 @ 1601

Last change on this file since 1601 was 1601, checked in by ctlod, 15 years ago

Doctor naming of OPA namelist variables , see ticket: #526

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 8.8 KB
Line 
1MODULE solpcg
2   !!======================================================================
3   !!                     ***  MODULE  solfet
4   !! Ocean solver :  preconditionned conjugate gradient solver
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   sol_pcg    : preconditionned conjugate gradient solver
9   !!----------------------------------------------------------------------
10   USE oce             ! ocean dynamics and tracers variables
11   USE dom_oce         ! ocean space and time domain variables
12   USE sol_oce         ! ocean solver variables
13   USE lib_mpp         ! distributed memory computing
14   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
15   USE in_out_manager  ! I/O manager
16
17   IMPLICIT NONE
18   PRIVATE
19
20   PUBLIC   sol_pcg    !
21
22   !! * Substitutions
23#  include "vectopt_loop_substitute.h90"
24   !!----------------------------------------------------------------------
25   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
26   !! $Id$
27   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
28   !!----------------------------------------------------------------------
29CONTAINS
30
31   SUBROUTINE sol_pcg( kindic )
32      !!----------------------------------------------------------------------
33      !!                  ***  ROUTINE sol_pcg  ***
34      !!                   
35      !! ** Purpose :   Solve the ellipic equation for the transport
36      !!      divergence system  using a diagonal preconditionned
37      !!      conjugate gradient method.
38      !!
39      !! ** Method  :   Diagonal preconditionned conjugate gradient method.
40      !!      the algorithm is multitasked. (case of 5 points matrix)
41      !!      define              pa  = q^-1 * a
42      !!                        pgcb  = q^-1 * gcb
43      !!                 < . ; . >_q  = ( . )^t q ( . )
44      !!      where q is the preconditioning matrix = diagonal matrix of the
45      !!                                              diagonal elements of a
46      !!      Initialization  :
47      !!         x(o) = gcx
48      !!         r(o) = d(o) = pgcb - pa.x(o)
49      !!         rr(o)= < r(o) , r(o) >_q
50      !!      Iteration 1     :
51      !!         standard PCG algorithm
52      !!      Iteration n > 1 :
53      !!         s(n)   = pa.r(n)
54      !!         gam(n) = < r(n) , r(n) >_q
55      !!         del(n) = < r(n) , s(n) >_q
56      !!         bet(n) = gam(n) / gam(n-1)
57      !!         d(n)   = r(n) + bet(n) d(n-1)
58      !!         z(n)   = s(n) + bet(n) z(n-1)
59      !!         sig(n) = del(n) - bet(n)*bet(n)*sig(n-1)
60      !!         alp(n) = gam(n) / sig(n)
61      !!         x(n+1) = x(n) + alp(n) d(n)
62      !!         r(n+1) = r(n) - alp(n) z(n)
63      !!      Convergence test :
64      !!         rr(n+1) / < gcb , gcb >_q   =< epsr
65      !!
66      !! ** Action : - niter  : solver number of iteration done
67      !!             - res    : solver residu reached
68      !!             - gcx()  : solution of the elliptic system
69      !!
70      !! References :
71      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
72      !!      D Azevedo et al. 1993, Computer Science Technical Report, Tennessee U.
73      !!
74      !! History :
75      !!        !  90-10  (G. Madec)  Original code
76      !!        !  91-11  (G. Madec)
77      !!        !  93-04  (M. Guyon)  loops and suppress pointers
78      !!        !  95-09  (M. Imbard, J. Escobar)  mpp exchange
79      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
80      !!        !  96-11  (A. Weaver)  correction to preconditioning
81      !!   8.5  !  02-08  (G. Madec)  F90: Free form
82      !!        !  08-01  (R. Benshila) mpp optimization
83      !!----------------------------------------------------------------------
84      INTEGER, INTENT( inout ) ::   kindic   ! solver indicator, < 0 if the conver-
85      !                                      ! gence is not reached: the model is
86      !                                      ! stopped in step
87      !                                      ! set to zero before the call of solpcg
88      !!
89      INTEGER ::   ji, jj, jn                ! dummy loop indices
90      REAL(wp) ::  zgcad                     ! temporary scalars
91      REAL(wp), DIMENSION(2) :: zsum
92      REAL(wp), DIMENSION(jpi,jpj) :: zgcr
93      !!----------------------------------------------------------------------
94
95      ! Initialization of the algorithm with standard PCG
96      ! -------------------------------------------------
97
98      CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
99
100      ! gcr   = gcb-a.gcx
101      ! gcdes = gcr
102
103      DO jj = 2, jpjm1
104         DO ji = fs_2, fs_jpim1   ! vector opt.
105            zgcad = bmask(ji,jj) * ( gcb(ji,jj  ) -                gcx(ji  ,jj  )   &
106               &                                  - gcp(ji,jj,1) * gcx(ji  ,jj-1)   &
107               &                                  - gcp(ji,jj,2) * gcx(ji-1,jj  )   &
108               &                                  - gcp(ji,jj,3) * gcx(ji+1,jj  )   &
109               &                                  - gcp(ji,jj,4) * gcx(ji  ,jj+1)   )
110            gcr  (ji,jj) = zgcad
111            gcdes(ji,jj) = zgcad
112         END DO
113      END DO
114
115      ! rnorme = (gcr,gcr)
116      rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
117      IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
118
119      CALL lbc_lnk( gcdes, c_solver_pt, 1. )   ! lateral boundary condition
120
121      ! gccd = matrix . gcdes
122      DO jj = 2, jpjm1
123         DO ji = fs_2, fs_jpim1   ! vector opt.
124            gccd(ji,jj) = bmask(ji,jj)*( gcdes(ji,jj)   &
125               &        +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
126               &        +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   )
127         END DO
128      END DO 
129
130      ! alph = (gcr,gcr)/(gcdes,gccd)
131      radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
132      IF( lk_mpp )   CALL mpp_sum( radd )   ! sum over the global domain
133      alph = rnorme /radd
134
135      ! gcx = gcx + alph * gcdes
136      ! gcr = gcr - alph * gccd
137      DO jj = 2, jpjm1
138         DO ji = fs_2, fs_jpim1   ! vector opt.
139            gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
140            gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
141         END DO
142      END DO
143
144      ! Algorithm wtih Eijkhout rearrangement
145      ! -------------------------------------
146       
147      !                                                !================
148      DO jn = 1, nn_nmax                               ! Iterative loop
149         !                                             !================
150
151         CALL lbc_lnk( gcr, c_solver_pt, 1. )   ! lateral boundary condition
152
153         ! zgcr = matrix . gcr
154         DO jj = 2, jpjm1
155            DO ji = fs_2, fs_jpim1   ! vector opt.
156               zgcr(ji,jj) = bmask(ji,jj)*( gcr(ji,jj)   &
157                  &        +gcp(ji,jj,1)*gcr(ji,jj-1)+gcp(ji,jj,2)*gcr(ji-1,jj)   &
158                  &        +gcp(ji,jj,4)*gcr(ji,jj+1)+gcp(ji,jj,3)*gcr(ji+1,jj)   )
159            END DO
160         END DO
161 
162         ! rnorme = (gcr,gcr)
163         rr = rnorme
164         zsum(1) = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
165
166         ! zgcad = (zgcr,gcr)
167         zsum(2) = SUM( gcr(2:jpim1,2:jpjm1) * gcdmat(2:jpim1,2:jpjm1) * zgcr(2:jpim1,2:jpjm1) )
168
169         IF( lk_mpp )   CALL mpp_sum( zsum, 2 )   ! sum over the global domain
170         rnorme = zsum(1) 
171         zgcad  = zsum(2)
172
173         ! test of convergence
174         IF( rnorme < epsr .OR. jn == nn_nmax ) THEN
175            res = SQRT( rnorme )
176            niter = jn
177            ncut = 999
178         ENDIF
179
180         ! beta = (rk+1,rk+1)/(rk,rk)
181         beta = rnorme / rr
182         radd = zgcad - beta*beta*radd
183         alph = rnorme / radd
184
185         ! gcx = gcx + alph * gcdes
186         ! gcr = gcr - alph * gccd
187         DO jj = 2, jpjm1
188            DO ji = fs_2, fs_jpim1   ! vector opt.
189               gcdes(ji,jj) = gcr (ji,jj) + beta * gcdes(ji,jj) 
190               gccd (ji,jj) = zgcr(ji,jj) + beta * gccd (ji,jj) 
191               gcx  (ji,jj) = gcx (ji,jj) + alph * gcdes(ji,jj) 
192               gcr  (ji,jj) = gcr (ji,jj) - alph * gccd (ji,jj) 
193            END DO
194         END DO
195       
196         ! indicator of non-convergence or explosion
197         IF( jn == nn_nmax .OR. SQRT(epsr)/eps > 1.e+20 ) kindic = -2
198         IF( ncut == 999 ) GOTO 999
199
200         !                                             !================
201      END DO                                           !    End Loop
202      !                                                !================
203     
204999   CONTINUE
205     
206     
207      ! Output in gcx with lateral b.c. applied
208      ! ---------------------------------------
209     
210      CALL lbc_lnk( gcx, c_solver_pt, 1. )
211     
212   END SUBROUTINE sol_pcg
213
214   !!=====================================================================
215END MODULE solpcg
Note: See TracBrowser for help on using the repository browser.