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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SOL – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90 @ 4460

Last change on this file since 4460 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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