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

source: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90 @ 3340

Last change on this file since 3340 was 3340, checked in by sga, 12 years ago

NEMO branch dev_r3337_NOCS10_ICB: add changes to ocean code to allow interface to iceberg code

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