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 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 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   !! * Modules used
11   USE oce             ! ocean dynamics and tracers variables
12   USE dom_oce         ! ocean space and time domain variables
13   USE sol_oce         ! ocean solver variables
14   USE lib_mpp         ! distributed memory computing
15   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
16
17   IMPLICIT NONE
18   PRIVATE
19
20   !! * Routine accessibility
21   PUBLIC sol_pcg              ! ???
22
23   !! * Substitutions
24#  include "vectopt_loop_substitute.h90"
25   !!----------------------------------------------------------------------
26
27CONTAINS
28
29   SUBROUTINE sol_pcg( kindic )
30      !!----------------------------------------------------------------------
31      !!                  ***  ROUTINE sol_pcg  ***
32      !!                   
33      !! ** Purpose :   Solve the ellipic equation for the barotropic stream
34      !!      function system (default option) or the transport divergence
35      !!      system ("key_dynspg_fsc") using a diagonal preconditionned
36      !!      conjugate gradient method.
37      !!      In the former case, the barotropic stream function trend has a
38      !!      zero boundary condition along all coastlines (i.e. continent
39      !!      as well as islands) while in the latter the boundary condition
40      !!      specification is not required.
41      !!
42      !! ** Method  :   Diagonal preconditionned conjugate gradient method.
43      !!      the algorithm is multitasked. (case of 5 points matrix)
44      !!      define              pa  = q^-1 * a
45      !!                        pgcb  = q^-1 * gcb
46      !!                 < . ; . >_q  = ( . )^t q ( . )
47      !!      where q is the preconditioning matrix = diagonal matrix of the
48      !!                                              diagonal elements of a
49      !!      Initialization:
50      !!         x(o) = gcx
51      !!         r(o) = d(o) = pgcb - pa.x(o)
52      !!         rr(o)= < r(o) , r(o) >_q
53      !!      Iteration n   :
54      !!         z(n)   = pa.d(n)
55      !!         alp(n) = rr(n) / < z(n) , d(n) >_q
56      !!         x(n+1) = x(n) + alp(n) d(n)
57      !!         r(n+1) = r(n) - alp(n) z(n)
58      !!         rr(n+1)= < r(n+1) , r(n+1) >_q
59      !!         bet(n) = rr(n+1) / rr(n)
60      !!         r(n+1) = r(n+1) + bet(n+1) d(n)
61      !!      Convergence test :
62      !!         rr(n+1) / < gcb , gcb >_q   =< epsr
63      !!
64      !! ** Action : - niter  : solver number of iteration done
65      !!             - res    : solver residu reached
66      !!             - gcx()  : solution of the elliptic system
67      !!
68      !! References :
69      !!      Madec et al. 1988, Ocean Modelling, issue 78, 1-6.
70      !!
71      !! History :
72      !!        !  90-10  (G. Madec)  Original code
73      !!        !  91-11  (G. Madec)
74      !!        !  93-04  (M. Guyon)  loops and suppress pointers
75      !!        !  95-09  (M. Imbard, J. Escobar)  mpp exchange
76      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
77      !!        !  96-11  (A. Weaver)  correction to preconditioning
78      !!   8.5  !  02-08  (G. Madec)  F90: Free form
79      !!----------------------------------------------------------------------
80      !! * Arguments
81      INTEGER, INTENT( inout ) ::   kindic   ! solver indicator, < 0 if the conver-
82      !                                      ! gence is not reached: the model is
83      !                                      ! stopped in step
84      !                                      ! set to zero before the call of solpcg
85
86      !! * Local declarations
87      INTEGER ::   ji, jj, jn                ! dummy loop indices
88      REAL(wp) ::   zgcad                    ! temporary scalars
89      !!----------------------------------------------------------------------
90
91      !                                                !================
92      DO jn = 1, nmax                                  ! Iterative loop
93         !                                             !================
94
95         !,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
96         
97         IF( jn == 1 ) THEN
98           
99            ! 1.0 Initialization of the algorithm
100            ! -----------------------------------
101           
102#if defined key_dynspg_fsc
103#   if defined key_mpp
104            ! Mpp: export boundary values to neighbouring processors
105            CALL lbc_lnk( gcx, 'S', 1. )
106#   else
107            !   mono- or macro-tasking: W-point, >0, 2D array, no slab
108            CALL lbc_lnk( gcx, 'T', 1. )
109#   endif
110#else
111#   if defined key_mpp
112            ! ... Mpp: export boundary values to neighbouring processors
113            CALL lbc_lnk( gcx, 'G', 1. )
114#   else
115            !   ... mono- or macro-tasking: F-point, >0, 2D array, no slab
116            CALL lbc_lnk( gcx, 'F', 1. )
117#   endif
118#endif
119
120            !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,,
121           
122            ! gcr   = gcb-a.gcx
123            ! gcdes = gsr
124           
125            DO jj = 2, jpjm1
126               DO ji = fs_2, fs_jpim1   ! vector opt.
127                  zgcad = bmask(ji,jj)*(                         &
128                    gcb(ji,jj  ) -              gcx(ji  ,jj  )   &
129                                 - gcp(ji,jj,1)*gcx(ji  ,jj-1)   &
130                                 - gcp(ji,jj,2)*gcx(ji-1,jj  )   &
131                                 - gcp(ji,jj,3)*gcx(ji+1,jj  )   &
132                                 - gcp(ji,jj,4)*gcx(ji  ,jj+1)   )
133                  gcr  (ji,jj) = zgcad
134                  gcdes(ji,jj) = zgcad
135               END DO
136            END DO
137           
138            !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,,
139           
140            rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
141
142#if defined key_mpp
143            ! Mpp: sum over all the global domain
144            CALL mpp_sum( rnorme )
145#endif
146            rr = rnorme
147
148        ENDIF
149        !,,,,,,,,,,,,,,,,,,,,,,,,synchro ,,,,,,,,,,,,,,,,,,,,,,,
150       
151       
152        ! 1.1 Algorithm
153        ! -------------
154       
155        ! boundary condition on gcdes (only cyclic bc are required)
156#if defined key_dynspg_fsc
157#   if defined key_mpp
158        !   Mpp: export boundary values to neighbouring processors
159        CALL lbc_lnk( gcdes, 'S', 1. )
160#   else
161        !   mono- or macro-tasking: W-point, >0, 2D array, no slab
162        CALL lbc_lnk( gcdes, 'T', 1. )
163#   endif
164#else
165#   if defined key_mpp
166        !   Mpp: export boundary values to neighbouring processors
167        CALL lbc_lnk( gcdes, 'G', 1. )
168#   else
169        !   mono- or macro-tasking: F-point, >0, 2D array, no slab
170        CALL lbc_lnk( gcdes, 'F', 1. )
171#   endif
172#endif
173       
174        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
175       
176        ! ... gccd = matrix . gcdes
177        DO jj = 2, jpjm1
178           DO ji = fs_2, fs_jpim1   ! vector opt.
179              gccd(ji,jj) = bmask(ji,jj)*(   &
180                 gcdes(ji,jj)   &
181                +gcp(ji,jj,1)*gcdes(ji,jj-1)+gcp(ji,jj,2)*gcdes(ji-1,jj)   &
182                +gcp(ji,jj,4)*gcdes(ji,jj+1)+gcp(ji,jj,3)*gcdes(ji+1,jj)   &
183                )
184           END DO
185        END DO
186
187        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
188       
189        ! alph = (gcr,gcr)/(gcdes,gccd)
190
191        radd = SUM(  gcdes(:,:) * gcdmat(:,:) * gccd(:,:)  )
192
193#if defined key_mpp
194        ! Mpp: sum over all the global domain
195        CALL mpp_sum( radd )
196#endif
197        alph = rr / radd
198       
199        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
200       
201        ! gcx = gcx + alph * gcdes
202        ! gcr = gcr - alph * gccd
203        DO jj = 2, jpjm1
204           DO ji = fs_2, fs_jpim1   ! vector opt.
205              gcx(ji,jj) = bmask(ji,jj) * ( gcx(ji,jj) + alph * gcdes(ji,jj) )
206              gcr(ji,jj) = bmask(ji,jj) * ( gcr(ji,jj) - alph * gccd (ji,jj) )
207           END DO
208        END DO
209       
210        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
211       
212        ! rnorme = (gcr,gcr)
213
214        rnorme = SUM(  gcr(:,:) * gcdmat(:,:) * gcr(:,:)  )
215
216#if defined key_mpp
217        ! Mpp: sum over all the global domain
218        CALL  mpp_sum( rnorme )
219#endif
220       
221        ! test of convergence
222        IF ( rnorme < epsr .OR. jn == nmax ) THEN
223           res = SQRT( rnorme )
224           niter = jn
225           ncut = 999
226        ENDIF
227       
228        ! beta = (rk+1,rk+1)/(rk,rk)
229        beta = rnorme / rr
230        rr   = rnorme
231
232        !,,,,,,,,,,,,,,,,,,,,,,,,synchro if macrotasking,,,,,,,,,,,,,,,,,,,,,,,
233
234        ! indicator of non-convergence or explosion
235        IF( jn == nmax .OR. sqrt(epsr)/eps > 1.e+20 ) kindic = -2
236        IF( ncut == 999 ) GOTO 999
237
238        ! gcdes = gcr + beta * gcdes
239        DO jj = 2, jpjm1
240           DO ji = fs_2, fs_jpim1   ! vector opt.
241              gcdes(ji,jj) = bmask(ji,jj)*( gcr(ji,jj) + beta * gcdes(ji,jj) )
242           END DO
243        END DO
244       
245        !                                             !================
246     END DO                                           !    End Loop
247     !                                                !================
248     
249999  CONTINUE
250     
251     
252     !  2. Output in gcx with lateral b.c. applied
253     !  ------------------------------------------
254     
255     ! boundary conditions   !!bug ???
256#if defined key_mpp
257     ! Mpp: export boundary values to neighbouring processors
258# if defined key_dynspg_fsc
259     CALL lbc_lnk( gcx, 'S', 1. )
260# else
261     CALL lbc_lnk( gcx, 'G', 1. )
262# endif
263#else
264     IF ( nperio /= 0 ) THEN
265# if defined key_dynspg_fsc
266        ! mono- or macro-tasking: W-point, >0, 2D array, no slab
267        CALL lbc_lnk( gcx, 'T', 1. )
268# else
269        ! mono- or macro-tasking: F-point, >0, 2D array, no slab
270        CALL lbc_lnk( gcx, 'F', 1. )
271# endif
272     ENDIF
273#endif
274     
275   END SUBROUTINE sol_pcg
276
277   !!=====================================================================
278END MODULE solpcg
Note: See TracBrowser for help on using the repository browser.