source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90 @ 4153

Last change on this file since 4153 was 4153, checked in by cetlod, 7 years ago

dev_LOCEAN_2013: merge in trunk changes between r3940 and r4028, see ticket #1169

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