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 in branches/devukmo2010/NEMO/OPA_SRC/SOL – NEMO

source: branches/devukmo2010/NEMO/OPA_SRC/SOL/solmat.F90 @ 2133

Last change on this file since 2133 was 2133, checked in by rfurner, 14 years ago

mistake in if statements in solmat corrected (issue was due to merge of bdy and trunk update from 3.2)

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