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

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

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