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.
solmat.F90 on Ticket #688 – Attachment – NEMO

Ticket #688: solmat.F90

File solmat.F90, 15.5 KB (added by molines, 14 years ago)
Line 
1MODULE solmat
2   !!======================================================================
3   !!                       ***  MODULE  solmat  ***
4   !! solver       : construction of the matrix
5   !!======================================================================
6   !! History :   1.0  ! 1988-04  (G. Madec)  Original code
7   !!                  ! 1993-03  (M. Guyon)  symetrical conditions
8   !!                  ! 1993-06  (M. Guyon)  suppress pointers
9   !!                  ! 1996-05  (G. Madec)  merge sor and pcg formulations
10   !!                  ! 1996-11  (A. Weaver)  correction to preconditioning
11   !!   NEMO      1.0  ! 1902-08  (G. Madec)  F90: Free form
12   !!              -   ! 1902-11  (C. Talandier, A-M. Treguier) Free surface & Open boundaries
13   !!             2.0  ! 2005-09  (R. Benshila)  add sol_exd for extra outer halo
14   !!              -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
15   !!             3.2  ! 2009-06  (S. Masson)  distributed restart using iom
16   !!              -   ! 2009-07  (R. Benshila)  suppression of rigid-lid option
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   sol_mat : Construction of the matrix of used by the elliptic solvers
21   !!   sol_exd :
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and active tracers
24   USE dom_oce         ! ocean space and time domain
25   USE sol_oce         ! ocean solver
26   USE phycst          ! physical constants
27   USE obc_oce         ! ocean open boundary conditions
28   USE lbclnk          ! lateral boudary conditions
29   USE lib_mpp         ! distributed memory computing
30   USE in_out_manager  ! I/O manager
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   sol_mat    ! routine called by inisol.F90
36
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
39   !! $Id: solmat.F90 1876 2010-05-18 15:35:35Z rblod $
40   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42
43CONTAINS
44
45   SUBROUTINE sol_mat( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE sol_mat  ***
48      !!
49      !! ** Purpose :   Construction of the matrix of used by the elliptic
50      !!              solvers (either sor or pcg methods).
51      !!
52      !! ** Method  :   The matrix is built for the divergence of the transport
53      !!              system. a diagonal preconditioning matrix is also defined.
54      !!
55      !! ** Action  : - gcp    : extra-diagonal elements of the matrix
56      !!              - gcdmat : preconditioning matrix (diagonal elements)
57      !!              - gcdprc : inverse of the preconditioning matrix
58      !!----------------------------------------------------------------------
59      INTEGER, INTENT(in) :: kt
60      !!
61      INTEGER ::   ji, jj                    ! dummy loop indices
62      REAL(wp) ::   zcoefs, zcoefw, zcoefe, zcoefn  ! temporary scalars
63      REAL(wp) ::   z2dt, zcoef
64      !!----------------------------------------------------------------------
65
66     
67      ! 1. Construction of the matrix
68      ! -----------------------------
69      zcoef = 0.e0                          ! initialize to zero
70      gcp(:,:,1) = 0.e0
71      gcp(:,:,2) = 0.e0
72      gcp(:,:,3) = 0.e0
73      gcp(:,:,4) = 0.e0
74      !
75      gcdprc(:,:) = 0.e0
76      gcdmat(:,:) = 0.e0
77      !
78      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt = rdt
79      ELSE                                        ;   z2dt = 2. * rdt
80      ENDIF
81
82#if defined key_dynspg_flt 
83#   if ! defined key_obc
84
85      DO jj = 2, jpjm1                      ! matrix of free surface elliptic system
86         DO ji = 2, jpim1
87            zcoef = z2dt * z2dt * grav * bmask(ji,jj)
88            zcoefs = -zcoef * hv(ji  ,jj-1) * e1v(ji  ,jj-1) / e2v(ji  ,jj-1)    ! south coefficient
89            zcoefw = -zcoef * hu(ji-1,jj  ) * e2u(ji-1,jj  ) / e1u(ji-1,jj  )    ! west coefficient
90            zcoefe = -zcoef * hu(ji  ,jj  ) * e2u(ji  ,jj  ) / e1u(ji  ,jj  )    ! east coefficient
91            zcoefn = -zcoef * hv(ji  ,jj  ) * e1v(ji  ,jj  ) / e2v(ji  ,jj  )    ! north coefficient
92            gcp(ji,jj,1) = zcoefs
93            gcp(ji,jj,2) = zcoefw
94            gcp(ji,jj,3) = zcoefe
95            gcp(ji,jj,4) = zcoefn
96            gcdmat(ji,jj) = e1t(ji,jj) * e2t(ji,jj) * bmask(ji,jj)    &          ! diagonal coefficient
97               &          - zcoefs -zcoefw -zcoefe -zcoefn
98         END DO
99      END DO
100#   else
101    IF ( Agrif_Root() ) THEN
102      DO jj = 2, jpjm1                      ! matrix of free surface elliptic system with open boundaries
103         DO ji = 2, jpim1
104            zcoef = z2dt * z2dt * grav * bmask(ji,jj)
105            !                                    ! south coefficient
106            IF( lp_obc_south .AND. ( jj == njs0p1 ) ) THEN
107               zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)*(1.-vsmsk(ji,1))
108            ELSE
109               zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)
110            END IF
111            gcp(ji,jj,1) = zcoefs
112            !
113            !                                    ! west coefficient
114            IF( lp_obc_west  .AND. ( ji == niw0p1 ) ) THEN
115               zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)*(1.-uwmsk(jj,1))
116            ELSE
117               zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)
118            END IF
119            gcp(ji,jj,2) = zcoefw
120            !
121            !                                    ! east coefficient
122            IF( lp_obc_east  .AND. ( ji == nie0 ) ) THEN
123               zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)*(1.-uemsk(jj,1))
124            ELSE
125               zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)
126            END IF
127            gcp(ji,jj,3) = zcoefe
128            !
129            !                                    ! north coefficient
130            IF( lp_obc_north .AND. ( jj == njn0 ) ) THEN
131               zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)*(1.-vnmsk(ji,1))
132            ELSE
133               zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)
134            END IF
135            gcp(ji,jj,4) = zcoefn
136            !
137            !                                    ! diagonal coefficient
138            gcdmat(ji,jj) = e1t(ji,jj)*e2t(ji,jj)*bmask(ji,jj)   &
139               &            - zcoefs -zcoefw -zcoefe -zcoefn
140         END DO
141      END DO
142    ELSE
143      DO jj = 2, jpjm1                      ! matrix of free surface elliptic system
144         DO ji = 2, jpim1
145            zcoef = z2dt * z2dt * grav * bmask(ji,jj)
146            zcoefs = -zcoef * hv(ji  ,jj-1) * e1v(ji  ,jj-1) / e2v(ji  ,jj-1)    ! south coefficient
147            zcoefw = -zcoef * hu(ji-1,jj  ) * e2u(ji-1,jj  ) / e1u(ji-1,jj  )    ! west coefficient
148            zcoefe = -zcoef * hu(ji  ,jj  ) * e2u(ji  ,jj  ) / e1u(ji  ,jj  )    ! east coefficient
149            zcoefn = -zcoef * hv(ji  ,jj  ) * e1v(ji  ,jj  ) / e2v(ji  ,jj  )    ! north coefficient
150            gcp(ji,jj,1) = zcoefs
151            gcp(ji,jj,2) = zcoefw
152            gcp(ji,jj,3) = zcoefe
153            gcp(ji,jj,4) = zcoefn
154            gcdmat(ji,jj) = e1t(ji,jj) * e2t(ji,jj) * bmask(ji,jj)    &          ! diagonal coefficient
155               &          - zcoefs -zcoefw -zcoefe -zcoefn
156         END DO
157      END DO
158    ENDIF
159#   endif
160#endif
161
162      IF( .NOT. Agrif_Root() ) THEN
163         !
164         IF( nbondi == -1 .OR. nbondi == 2 )   bmask(2     ,:     ) = 0.e0
165         IF( nbondi ==  1 .OR. nbondi == 2 )   bmask(nlci-1,:     ) = 0.e0
166         IF( nbondj == -1 .OR. nbondj == 2 )   bmask(:     ,2     ) = 0.e0
167         IF( nbondj ==  1 .OR. nbondj == 2 )   bmask(:     ,nlcj-1) = 0.e0
168         !
169         DO jj = 2, jpjm1
170            DO ji = 2, jpim1
171               zcoef = z2dt * z2dt * grav * bmask(ji,jj)
172               !  south coefficient
173               IF( ( nbondj == -1 .OR. nbondj == 2 ) .AND. ( jj == 3 ) ) THEN
174                  zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)*(1.-vmask(ji,jj-1,1))
175               ELSE
176                  zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)
177               END IF
178               gcp(ji,jj,1) = zcoefs
179               !
180               !  west coefficient
181               IF( ( nbondi == -1 .OR. nbondi == 2 ) .AND. ( ji == 3 )  ) THEN
182                  zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)*(1.-umask(ji-1,jj,1))
183               ELSE
184                  zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)
185               END IF
186               gcp(ji,jj,2) = zcoefw
187               !
188               !   east coefficient
189               IF( ( nbondi == 1 .OR. nbondi == 2 ) .AND. ( ji == nlci-2 ) ) THEN
190                  zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)*(1.-umask(ji,jj,1))
191               ELSE
192                  zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)
193               END IF
194               gcp(ji,jj,3) = zcoefe
195               !
196               !   north coefficient
197               IF( ( nbondj == 1 .OR. nbondj == 2 ) .AND. ( jj == nlcj-2 ) ) THEN
198                  zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)*(1.-vmask(ji,jj,1))
199               ELSE
200                  zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)
201               END IF
202               gcp(ji,jj,4) = zcoefn
203               !
204               ! diagonal coefficient
205               gcdmat(ji,jj) = e1t(ji,jj)*e2t(ji,jj)*bmask(ji,jj)   &
206                  &            - zcoefs -zcoefw -zcoefe -zcoefn
207            END DO
208         END DO
209         !
210      ENDIF
211
212      ! 2. Boundary conditions
213      ! ----------------------
214     
215      ! Cyclic east-west boundary conditions
216      !     ji=2 is the column east of ji=jpim1 and reciprocally,
217      !     ji=jpim1 is the column west of ji=2
218      !     all the coef are already set to zero as bmask is initialized to
219      !     zero for ji=1 and ji=jpj in dommsk.
220     
221      ! Symetrical conditions
222      ! free surface: no specific action
223      ! bsf system: n-s gradient of bsf = 0 along j=2 (perhaps a bug !!!!!!)
224      ! the diagonal coefficient of the southern grid points must be modify to
225      ! account for the existence of the south symmetric bassin.
226     
227      ! North fold boundary condition
228      !     all the coef are already set to zero as bmask is initialized to
229      !     zero on duplicated lignes and portion of lignes
230     
231      ! 3. Preconditioned matrix
232      ! ------------------------
233     
234      ! SOR and PCG solvers
235      DO jj = 1, jpj
236         DO ji = 1, jpi
237            IF( bmask(ji,jj) /= 0.e0 )   gcdprc(ji,jj) = 1.e0 / gcdmat(ji,jj)
238         END DO
239      END DO
240         
241      gcp(:,:,1) = gcp(:,:,1) * gcdprc(:,:)
242      gcp(:,:,2) = gcp(:,:,2) * gcdprc(:,:)
243      gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:)
244      gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:)
245      IF( nn_solv == 2 )  gccd(:,:) = rn_sor * gcp(:,:,2)
246
247      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0) THEN
248         CALL lbc_lnk_e( gcp   (:,:,1), c_solver_pt, 1. )   ! lateral boundary conditions
249         CALL lbc_lnk_e( gcp   (:,:,2), c_solver_pt, 1. )   ! lateral boundary conditions
250         CALL lbc_lnk_e( gcp   (:,:,3), c_solver_pt, 1. )   ! lateral boundary conditions
251         CALL lbc_lnk_e( gcp   (:,:,4), c_solver_pt, 1. )   ! lateral boundary conditions
252         CALL lbc_lnk_e( gcdprc(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions
253         CALL lbc_lnk_e( gcdmat(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions         
254         IF( npolj /= 0 ) CALL sol_exd( gcp , c_solver_pt ) ! switch northernelements
255      END IF
256   
257      ! 4. Initialization the arrays used in pcg
258      ! ----------------------------------------
259      gcb  (:,:) = 0.e0
260      gcr  (:,:) = 0.e0
261      gcdes(:,:) = 0.e0
262      gccd (:,:) = 0.e0
263      !
264   END SUBROUTINE sol_mat
265
266
267   SUBROUTINE sol_exd( pt3d, cd_type )
268      !!----------------------------------------------------------------------
269      !!                  ***  routine sol_exd  ***
270      !!                 
271      !! ** Purpose :   Reorder gcb coefficient on the extra outer  halo
272      !!                at north fold in case of T or F pivot
273      !!
274      !! ** Method  :   Perform a circular permutation of the coefficients on
275      !!                the total area strictly above the pivot point,
276      !!                and on the semi-row of the pivot point   
277      !!----------------------------------------------------------------------
278      CHARACTER(len=1) , INTENT( in ) ::   cd_type   ! define the nature of pt2d array grid-points
279         !                                           !  = T , U , V , F , W
280         !                                           !  = S : T-point, north fold treatment
281         !                                           !  = G : F-point, north fold treatment
282         !                                           !  = I : sea-ice velocity at F-point with index shift
283      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4), INTENT(inout) ::   pt3d   ! 2D field to be treated
284      !!
285      INTEGER  ::   ji, jk   ! dummy loop indices
286      INTEGER  ::   iloc     ! temporary integers
287      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4) ::   ztab   ! 2D workspace
288      !!----------------------------------------------------------------------
289
290      ztab = pt3d
291
292      SELECT CASE ( npolj )            ! north fold type
293      !
294      CASE ( 3 , 4 )                        !==  T pivot  ==!
295         iloc = jpiglo/2 +1 
296         !   
297         SELECT CASE ( cd_type )
298         !
299         CASE ( 'T', 'S', 'U', 'W' )
300            DO jk = 1, 4
301               DO ji = 1-jpr2di, nlci+jpr2di
302                  pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
303               END DO
304            END DO
305            DO jk =1, 4
306               DO ji = nlci+jpr2di, 1-jpr2di,  -1
307                  IF( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) &
308                     & .OR. ( mi0(iloc) == jpi+1 ) ) EXIT
309                     pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
310               END DO
311            END DO
312            !
313         CASE ( 'F' ,'G' , 'I', 'V' )
314            DO jk =1, 4
315               DO ji = 1-jpr2di, nlci+jpr2di
316                  pt3d(ji,nlcj-1:nlcj+jpr2dj,jk) = ztab(ji,nlcj-1:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
317               END DO
318            END DO
319            !
320         END SELECT   ! cd_type
321          !
322      CASE ( 5 , 6 )                        !==  F pivot  ==!
323         iloc=jpiglo/2
324         !
325         SELECT CASE (cd_type )
326         !
327         CASE ( 'T'  ,'S', 'U', 'W')
328            DO jk =1, 4
329               DO ji = 1-jpr2di, nlci+jpr2di
330                  pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
331               END DO
332            END DO
333            !
334         CASE ( 'F' ,'G' , 'I', 'V' )
335            DO jk =1, 4
336               DO ji = 1-jpr2di, nlci+jpr2di
337                  pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
338               END DO
339            END DO
340            DO jk =1, 4
341               DO ji = nlci+jpr2di, 1-jpr2di,  -1
342                  IF( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) .OR. ( mi0(iloc) == jpi+1 ) )   EXIT
343                    pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
344               END DO
345            END DO
346            !
347         END SELECT   ! cd_type
348         !
349      END SELECT   ! npolj
350      !   
351   END SUBROUTINE sol_exd
352
353   !!======================================================================
354END MODULE solmat