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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/SOL – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90 @ 6137

Last change on this file since 6137 was 5870, checked in by acc, 9 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 14.3 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 bdy_oce         ! unstructured open boundary conditions
29   USE lbclnk          ! lateral boudary conditions
30   USE lib_mpp         ! distributed memory computing
31   USE c1d               ! 1D vertical configuration
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
88      DO jj = 2, jpjm1                      ! matrix of free surface elliptic system
89         DO ji = 2, jpim1
90            zcoef = z2dt * z2dt * grav * bmask(ji,jj)
91            zcoefs = -zcoef * hv(ji  ,jj-1) * e1_e2v(ji  ,jj-1)    ! south coefficient
92            zcoefw = -zcoef * hu(ji-1,jj  ) * e2_e1u(ji-1,jj  )    ! west coefficient
93            zcoefe = -zcoef * hu(ji  ,jj  ) * e2_e1u(ji  ,jj  )    ! east coefficient
94            zcoefn = -zcoef * hv(ji  ,jj  ) * e1_e2v(ji  ,jj  )    ! north coefficient
95            gcp(ji,jj,1) = zcoefs
96            gcp(ji,jj,2) = zcoefw
97            gcp(ji,jj,3) = zcoefe
98            gcp(ji,jj,4) = zcoefn
99            gcdmat(ji,jj) = e1e2t(ji,jj) * bmask(ji,jj)    &          ! diagonal coefficient
100               &          - zcoefs -zcoefw -zcoefe -zcoefn
101         END DO
102      END DO
103
104#  elif defined key_dynspg_flt && defined key_bdy 
105
106      !   defined gcdmat in the case of unstructured open boundaries
107      DO jj = 2, jpjm1
108         DO ji = 2, jpim1
109            zcoef = z2dt * z2dt * grav * bmask(ji,jj)
110
111            !  south coefficient
112            zcoefs = -zcoef * hv(ji,jj-1) * e1_e2v(ji,jj-1)
113            zcoefs = zcoefs * bdyvmask(ji,jj-1)
114            gcp(ji,jj,1) = zcoefs
115
116            !  west coefficient
117            zcoefw = -zcoef * hu(ji-1,jj) * e2_e1u(ji-1,jj)
118            zcoefw = zcoefw * bdyumask(ji-1,jj)
119            gcp(ji,jj,2) = zcoefw
120
121            !  east coefficient
122            zcoefe = -zcoef * hu(ji,jj) * e2_e1u(ji,jj)
123            zcoefe = zcoefe * bdyumask(ji,jj)
124            gcp(ji,jj,3) = zcoefe
125
126            !  north coefficient
127            zcoefn = -zcoef * hv(ji,jj) * e1_e2v(ji,jj)
128            zcoefn = zcoefn * bdyvmask(ji,jj)
129            gcp(ji,jj,4) = zcoefn
130
131            ! diagonal coefficient
132            gcdmat(ji,jj) = e1e2t(ji,jj)*bmask(ji,jj) - zcoefs -zcoefw -zcoefe -zcoefn
133         END DO
134      END DO
135
136#endif
137
138      IF( .NOT. Agrif_Root() ) THEN
139         !
140         IF( nbondi == -1 .OR. nbondi == 2 )   bmask(2     ,:     ) = 0.e0
141         IF( nbondi ==  1 .OR. nbondi == 2 )   bmask(nlci-1,:     ) = 0.e0
142         IF( nbondj == -1 .OR. nbondj == 2 )   bmask(:     ,2     ) = 0.e0
143         IF( nbondj ==  1 .OR. nbondj == 2 )   bmask(:     ,nlcj-1) = 0.e0
144         !
145         DO jj = 2, jpjm1
146            DO ji = 2, jpim1
147               zcoef = z2dt * z2dt * grav * bmask(ji,jj)
148               !  south coefficient
149               IF( ( nbondj == -1 .OR. nbondj == 2 ) .AND. ( jj == 3 ) ) THEN
150                  zcoefs = -zcoef * hv(ji,jj-1) * e1_e2v(ji,jj-1)*(1.-vmask(ji,jj-1,1))
151               ELSE
152                  zcoefs = -zcoef * hv(ji,jj-1) * e1_e2v(ji,jj-1)
153               END IF
154               gcp(ji,jj,1) = zcoefs
155               !
156               !  west coefficient
157               IF( ( nbondi == -1 .OR. nbondi == 2 ) .AND. ( ji == 3 )  ) THEN
158                  zcoefw = -zcoef * hu(ji-1,jj) * e2_e1u(ji-1,jj)*(1.-umask(ji-1,jj,1))
159               ELSE
160                  zcoefw = -zcoef * hu(ji-1,jj) * e2_e1u(ji-1,jj)
161               END IF
162               gcp(ji,jj,2) = zcoefw
163               !
164               !   east coefficient
165               IF( ( nbondi == 1 .OR. nbondi == 2 ) .AND. ( ji == nlci-2 ) ) THEN
166                  zcoefe = -zcoef * hu(ji,jj) * e2_e1u(ji,jj)*(1.-umask(ji,jj,1))
167               ELSE
168                  zcoefe = -zcoef * hu(ji,jj) * e2_e1u(ji,jj)
169               END IF
170               gcp(ji,jj,3) = zcoefe
171               !
172               !   north coefficient
173               IF( ( nbondj == 1 .OR. nbondj == 2 ) .AND. ( jj == nlcj-2 ) ) THEN
174                  zcoefn = -zcoef * hv(ji,jj) * e1_e2v(ji,jj)*(1.-vmask(ji,jj,1))
175               ELSE
176                  zcoefn = -zcoef * hv(ji,jj) * e1_e2v(ji,jj)
177               END IF
178               gcp(ji,jj,4) = zcoefn
179               !
180               ! diagonal coefficient
181               gcdmat(ji,jj) = e1e2t(ji,jj)*bmask(ji,jj) - zcoefs -zcoefw -zcoefe -zcoefn
182            END DO
183         END DO
184         !
185      ENDIF
186
187      ! 2. Boundary conditions
188      ! ----------------------
189     
190      ! Cyclic east-west boundary conditions
191      !     ji=2 is the column east of ji=jpim1 and reciprocally,
192      !     ji=jpim1 is the column west of ji=2
193      !     all the coef are already set to zero as bmask is initialized to
194      !     zero for ji=1 and ji=jpj in dommsk.
195     
196      ! Symetrical conditions
197      ! free surface: no specific action
198      ! bsf system: n-s gradient of bsf = 0 along j=2 (perhaps a bug !!!!!!)
199      ! the diagonal coefficient of the southern grid points must be modify to
200      ! account for the existence of the south symmetric bassin.
201     
202      ! North fold boundary condition
203      !     all the coef are already set to zero as bmask is initialized to
204      !     zero on duplicated lignes and portion of lignes
205     
206      ! 3. Preconditioned matrix
207      ! ------------------------
208     
209      ! SOR and PCG solvers
210      IF( lk_c1d ) CALL lbc_lnk( gcdmat, 'T', 1._wp ) ! 1D case bmask =/0  but gcdmat not define everywhere
211      DO jj = 1, jpj
212         DO ji = 1, jpi
213            IF( bmask(ji,jj) /= 0.e0 )   gcdprc(ji,jj) = 1.e0 / gcdmat(ji,jj)
214         END DO
215      END DO
216         
217      gcp(:,:,1) = gcp(:,:,1) * gcdprc(:,:)
218      gcp(:,:,2) = gcp(:,:,2) * gcdprc(:,:)
219      gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:)
220      gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:)
221      IF( nn_solv == 2 )  gccd(:,:) = rn_sor * gcp(:,:,2)
222
223      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0) THEN
224         CALL lbc_lnk_e( gcp   (:,:,1), c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
225         CALL lbc_lnk_e( gcp   (:,:,2), c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
226         CALL lbc_lnk_e( gcp   (:,:,3), c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
227         CALL lbc_lnk_e( gcp   (:,:,4), c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
228         CALL lbc_lnk_e( gcdprc(:,:)  , c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions
229         CALL lbc_lnk_e( gcdmat(:,:)  , c_solver_pt, 1., jpr2di, jpr2dj )   ! lateral boundary conditions         
230         IF( npolj /= 0 ) CALL sol_exd( gcp , c_solver_pt ) ! switch northernelements
231      END IF
232   
233      ! 4. Initialization the arrays used in pcg
234      ! ----------------------------------------
235      gcb  (:,:) = 0.e0
236      gcr  (:,:) = 0.e0
237      gcdes(:,:) = 0.e0
238      gccd (:,:) = 0.e0
239      !
240      IF( nn_timing == 1 )  CALL timing_stop('sol_mat')
241      !
242   END SUBROUTINE sol_mat
243
244
245   SUBROUTINE sol_exd( pt3d, cd_type )
246      !!----------------------------------------------------------------------
247      !!                  ***  routine sol_exd  ***
248      !!                 
249      !! ** Purpose :   Reorder gcb coefficient on the extra outer  halo
250      !!                at north fold in case of T or F pivot
251      !!
252      !! ** Method  :   Perform a circular permutation of the coefficients on
253      !!                the total area strictly above the pivot point,
254      !!                and on the semi-row of the pivot point   
255      !!----------------------------------------------------------------------
256      CHARACTER(len=1) , INTENT( in ) ::   cd_type   ! define the nature of pt2d array grid-points
257         !                                           !  = T , U , V , F , W
258         !                                           !  = S : T-point, north fold treatment
259         !                                           !  = G : F-point, north fold treatment
260         !                                           !  = I : sea-ice velocity at F-point with index shift
261      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4), INTENT(inout) ::   pt3d   ! 2D field to be treated
262      !!
263      INTEGER  ::   ji, jk   ! dummy loop indices
264      INTEGER  ::   iloc     ! local integers
265      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ztab   ! workspace allocated one for all
266      !!----------------------------------------------------------------------
267
268      IF( .NOT. ALLOCATED( ztab ) ) THEN
269         ALLOCATE( ztab(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4), STAT=iloc )
270         IF( lk_mpp    )   CALL mpp_sum ( iloc )
271         IF( iloc /= 0 )   CALL ctl_stop('STOP', 'sol_exd: failed to allocate array')
272      ENDIF
273     
274      ztab = pt3d
275
276      SELECT CASE ( npolj )            ! north fold type
277      !
278      CASE ( 3 , 4 )                        !==  T pivot  ==!
279         iloc = jpiglo/2 +1 
280         !   
281         SELECT CASE ( cd_type )
282         !
283         CASE ( 'T' , 'U', 'W' )
284            DO jk = 1, 4
285               DO ji = 1-jpr2di, nlci+jpr2di
286                  pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
287               END DO
288            END DO
289            DO jk =1, 4
290               DO ji = nlci+jpr2di, 1-jpr2di,  -1
291                  IF( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) &
292                     & .OR. ( mi0(iloc) == jpi+1 ) ) EXIT
293                     pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
294               END DO
295            END DO
296            !
297         CASE ( 'F' , 'I', 'V' )
298            DO jk =1, 4
299               DO ji = 1-jpr2di, nlci+jpr2di
300                  pt3d(ji,nlcj-1:nlcj+jpr2dj,jk) = ztab(ji,nlcj-1:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
301               END DO
302            END DO
303            !
304         END SELECT   ! cd_type
305          !
306      CASE ( 5 , 6 )                        !==  F pivot  ==!
307         iloc=jpiglo/2
308         !
309         SELECT CASE (cd_type )
310         !
311         CASE ( 'T' , 'U', 'W')
312            DO jk =1, 4
313               DO ji = 1-jpr2di, nlci+jpr2di
314                  pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
315               END DO
316            END DO
317            !
318         CASE ( 'F' , 'I', 'V' )
319            DO jk =1, 4
320               DO ji = 1-jpr2di, nlci+jpr2di
321                  pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
322               END DO
323            END DO
324            DO jk =1, 4
325               DO ji = nlci+jpr2di, 1-jpr2di,  -1
326                  IF( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) .OR. ( mi0(iloc) == jpi+1 ) )   EXIT
327                    pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
328               END DO
329            END DO
330            !
331         END SELECT   ! cd_type
332         !
333      END SELECT   ! npolj
334      !   
335   END SUBROUTINE sol_exd
336
337   !!======================================================================
338END MODULE solmat
Note: See TracBrowser for help on using the repository browser.