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 trunk/NEMO/OPA_SRC/SOL – NEMO

source: trunk/NEMO/OPA_SRC/SOL/solmat.F90 @ 896

Last change on this file since 896 was 784, checked in by rblod, 16 years ago

merge solsor and solsor_e, see ticket #45
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.4 KB
Line 
1MODULE solmat
2   !!======================================================================
3   !!                       ***  MODULE  solmat  ***
4   !! solver       : construction of the matrix
5   !!======================================================================
6   !! History :   1.0  !  88-04  (G. Madec)  Original code
7   !!                  !  93-03  (M. Guyon)  symetrical conditions
8   !!                  !  93-06  (M. Guyon)  suppress pointers
9   !!                  !  96-05  (G. Madec)  merge sor and pcg formulations
10   !!                  !  96-11  (A. Weaver)  correction to preconditioning
11   !!                  !  98-02  (M. Guyon)  FETI method
12   !!             8.5  !  02-08  (G. Madec)  F90: Free form
13   !!                  !  02-11  (C. Talandier, A-M. Treguier) Free surface & Open boundaries
14   !!             9.0  !  05-09  (R. Benshila)  add sol_exd for extra outer halo
15   !!             9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
16   !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
17   !!----------------------------------------------------------------------
18
19   !!----------------------------------------------------------------------
20   !!   sol_mat       : Construction of the matrix of used by the elliptic solvers
21   !!   fetsch        :
22   !!   fetmat        :
23   !!   fetstr        :
24   !!----------------------------------------------------------------------
25   !! * Modules used
26   USE oce             ! ocean dynamics and active tracers
27   USE dom_oce         ! ocean space and time domain
28   USE sol_oce         ! ocean solver
29   USE phycst          ! physical constants
30   USE obc_oce         ! ocean open boundary conditions
31   USE lbclnk          ! lateral boudary conditions
32   USE lib_mpp         ! distributed memory computing
33   USE in_out_manager  ! I/O manager
34
35   IMPLICIT NONE
36   PRIVATE
37
38   !! * Routine accessibility
39   PUBLIC sol_mat     ! routine called by inisol.F90
40   !!----------------------------------------------------------------------
41   !!   OPA 9.0 , LOCEAN-IPSL (2005)
42   !! $Header$
43   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
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, pcg or feti methods).
54      !!
55      !! ** Method  :   The matrix depends on the type of free surface:
56      !!       * lk_dynspg_rl=T: rigid lid formulation
57      !!      The matrix is built for the barotropic stream function system.
58      !!      a diagonal preconditioning matrix is also defined.
59      !!       * lk_dynspg_flt=T: free surface formulation
60      !!      The matrix is built for the divergence of the transport system
61      !!      a diagonal preconditioning matrix is also defined.
62      !!        Note that for feti solver (nsolv=3) a specific initialization
63      !!      is required (call to fetstr.F) for memory allocation and inter-
64      !!      face definition.
65      !!
66      !! ** Action  : - gcp    : extra-diagonal elements of the matrix
67      !!              - gcdmat : preconditioning matrix (diagonal elements)
68      !!              - gcdprc : inverse of the preconditioning matrix
69      !!----------------------------------------------------------------------
70      !! * Arguments
71      INTEGER, INTENT(in) :: kt
72
73      !! * Local declarations
74      INTEGER ::   ji, jj                    ! dummy loop indices
75      INTEGER ::   ii, ij, iiend, ijend      ! temporary integers
76      REAL(wp) ::   zcoefs, zcoefw, zcoefe, zcoefn  ! temporary scalars
77      REAL(wp) ::   z2dt, zcoef
78      !!----------------------------------------------------------------------
79
80      ! FETI method ( nsolv = 3)
81      ! memory allocation and interface definition for the solver
82
83      IF( nsolv == 3 )   CALL fetstr
84
85     
86      ! 1. Construction of the matrix
87      ! -----------------------------
88     
89      ! initialize to zero
90      zcoef = 0.e0
91      gcp(:,:,1) = 0.e0
92      gcp(:,:,2) = 0.e0
93      gcp(:,:,3) = 0.e0
94      gcp(:,:,4) = 0.e0
95     
96      gcdprc(:,:) = 0.e0
97      gcdmat(:,:) = 0.e0
98     
99      IF( neuler == 0 .AND. kt == nit000 ) THEN
100         z2dt = rdt
101      ELSE
102         z2dt = 2. * rdt
103      ENDIF
104
105#if defined key_dynspg_flt && ! defined key_obc
106!!cr      IF( lk_dynspg_flt .AND. .NOT.lk_obc ) THEN   !bug missing lk_dynspg_flt_atsk
107
108      ! defined the coefficients for free surface elliptic system
109
110      DO jj = 2, jpjm1
111         DO ji = 2, jpim1
112            zcoef = z2dt * z2dt * grav * rnu * bmask(ji,jj)
113            zcoefs = -zcoef * hv(ji  ,jj-1) * e1v(ji  ,jj-1) / e2v(ji  ,jj-1)    ! south coefficient
114            zcoefw = -zcoef * hu(ji-1,jj  ) * e2u(ji-1,jj  ) / e1u(ji-1,jj  )    ! west coefficient
115            zcoefe = -zcoef * hu(ji  ,jj  ) * e2u(ji  ,jj  ) / e1u(ji  ,jj  )    ! east coefficient
116            zcoefn = -zcoef * hv(ji  ,jj  ) * e1v(ji  ,jj  ) / e2v(ji  ,jj  )    ! north coefficient
117            gcp(ji,jj,1) = zcoefs
118            gcp(ji,jj,2) = zcoefw
119            gcp(ji,jj,3) = zcoefe
120            gcp(ji,jj,4) = zcoefn
121            gcdmat(ji,jj) = e1t(ji,jj) * e2t(ji,jj) * bmask(ji,jj)    &          ! diagonal coefficient
122               &          - zcoefs -zcoefw -zcoefe -zcoefn
123         END DO
124      END DO
125     
126#  elif defined key_dynspg_flt && defined key_obc
127!!cr      ELSEIF( lk_dynspg_flt .AND. lk_obc ) THEN     !bug missing lk_dynspg_flt_atsk
128
129      !   defined gcdmat in the case of open boundaries
130
131      DO jj = 2, jpjm1
132         DO ji = 2, jpim1
133            zcoef = z2dt * z2dt * grav * rnu * bmask(ji,jj)
134            !  south coefficient
135            IF( lp_obc_south .AND. ( jj == njs0p1 ) ) THEN
136               zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)*(1.-vsmsk(ji,1))
137            ELSE
138               zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)
139            END IF
140            gcp(ji,jj,1) = zcoefs
141
142            !  west coefficient
143            IF( lp_obc_west  .AND. ( ji == niw0p1 ) ) THEN
144               zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)*(1.-uwmsk(jj,1))
145            ELSE
146               zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)
147            END IF
148            gcp(ji,jj,2) = zcoefw
149
150            !   east coefficient
151            IF( lp_obc_east  .AND. ( ji == nie0 ) ) THEN
152               zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)*(1.-uemsk(jj,1))
153            ELSE
154               zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)
155            END IF
156            gcp(ji,jj,3) = zcoefe
157
158            !   north coefficient
159            IF( lp_obc_north .AND. ( jj == njn0 ) ) THEN
160               zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)*(1.-vnmsk(ji,1))
161            ELSE
162               zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)
163            END IF
164            gcp(ji,jj,4) = zcoefn
165
166            ! diagonal coefficient
167            gcdmat(ji,jj) = e1t(ji,jj)*e2t(ji,jj)*bmask(ji,jj) &
168                            - zcoefs -zcoefw -zcoefe -zcoefn
169         END DO
170      END DO
171
172#  else
173!!cr      ELSE
174
175      !   defined the coefficients for bsf elliptic system
176     
177      DO jj = 2, jpjm1
178         DO ji = 2, jpim1
179            zcoefs = -hur(ji  ,jj  ) * e1u(ji  ,jj  ) / e2u(ji  ,jj  ) * bmask(ji,jj)   ! south coefficient
180            zcoefw = -hvr(ji  ,jj  ) * e2v(ji  ,jj  ) / e1v(ji  ,jj  ) * bmask(ji,jj)   ! west coefficient
181            zcoefe = -hvr(ji+1,jj  ) * e2v(ji+1,jj  ) / e1v(ji+1,jj  ) * bmask(ji,jj)   ! east coefficient
182            zcoefn = -hur(ji  ,jj+1) * e1u(ji  ,jj+1) / e2u(ji  ,jj+1) * bmask(ji,jj)   ! north coefficient
183            gcp(ji,jj,1) = zcoefs
184            gcp(ji,jj,2) = zcoefw
185            gcp(ji,jj,3) = zcoefe
186            gcp(ji,jj,4) = zcoefn
187            gcdmat(ji,jj) = -zcoefs -zcoefw -zcoefe -zcoefn                             ! diagonal coefficient
188         END DO
189      END DO
190     
191!!cr  ENDIF
192#endif
193#if defined key_agrif
194       IF (.NOT.AGRIF_ROOT()) THEN
195       
196       IF ( (nbondi == -1)  .OR. (nbondi == 2) ) bmask(2,:)=0.
197       IF ( (nbondi ==  1)  .OR. (nbondi == 2) ) bmask(nlci-1,:)=0.
198       IF ( (nbondj == -1)  .OR. (nbondj == 2) ) bmask(:,2)=0.
199       IF ( (nbondj ==  1)  .OR. (nbondj == 2) ) bmask(:,nlcj-1)=0.
200
201      DO jj = 2, jpjm1
202         DO ji = 2, jpim1
203            zcoef = z2dt * z2dt * grav * rnu * bmask(ji,jj)
204            !  south coefficient
205            IF( ((nbondj == -1)  .OR. (nbondj == 2)) .AND. ( jj == 3 ) ) THEN
206               zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)*(1.-vmask(ji,jj-1,1))
207            ELSE
208               zcoefs = -zcoef * hv(ji,jj-1) * e1v(ji,jj-1)/e2v(ji,jj-1)
209            END IF
210            gcp(ji,jj,1) = zcoefs
211
212            !  west coefficient
213       IF( ( (nbondi == -1)  .OR. (nbondi == 2) ) .AND. ( ji == 3 )  ) THEN
214               zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)*(1.-umask(ji-1,jj,1))
215            ELSE
216               zcoefw = -zcoef * hu(ji-1,jj) * e2u(ji-1,jj)/e1u(ji-1,jj)
217            END IF
218            gcp(ji,jj,2) = zcoefw
219
220            !   east coefficient
221            IF( ((nbondi == 1)  .OR. (nbondi == 2)) .AND. ( ji == nlci-2 ) ) THEN
222               zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)*(1.-umask(ji,jj,1))
223            ELSE
224               zcoefe = -zcoef * hu(ji,jj) * e2u(ji,jj)/e1u(ji,jj)
225            END IF
226            gcp(ji,jj,3) = zcoefe
227
228            !   north coefficient
229            IF( ((nbondj == 1)  .OR. (nbondj == 2)) .AND. ( jj == nlcj-2 ) ) THEN
230               zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)*(1.-vmask(ji,jj,1))
231            ELSE
232               zcoefn = -zcoef * hv(ji,jj) * e1v(ji,jj)/e2v(ji,jj)
233            END IF
234            gcp(ji,jj,4) = zcoefn
235
236            ! diagonal coefficient
237            gcdmat(ji,jj) = e1t(ji,jj)*e2t(ji,jj)*bmask(ji,jj) &
238                            - zcoefs -zcoefw -zcoefe -zcoefn
239         END DO
240      END DO
241     
242       ENDIF
243#endif
244
245      ! 2. Boundary conditions
246      ! ----------------------
247     
248      ! Cyclic east-west boundary conditions
249      !     ji=2 is the column east of ji=jpim1 and reciprocally,
250      !     ji=jpim1 is the column west of ji=2
251      !     all the coef are already set to zero as bmask is initialized to
252      !     zero for ji=1 and ji=jpj in dommsk.
253     
254      ! Symetrical conditions
255      ! free surface: no specific action
256      ! bsf system: n-s gradient of bsf = 0 along j=2 (perhaps a bug !!!!!!)
257      ! the diagonal coefficient of the southern grid points must be modify to
258      ! account for the existence of the south symmetric bassin.
259     
260!!cr      IF( .NOT.lk_dynspg_flt ) THEN   !bug missing lk_dynspg_flt_atsk
261#if ! defined key_dynspg_flt
262      IF( nperio == 2 ) THEN
263         DO ji = 1, jpi
264            IF( bmask(ji,2) /= 0.e0 ) THEN
265               zcoefs = - hur(ji,2)*e1u(ji,2)/e2u(ji,2)
266               gcdmat(ji,2) = gcdmat(ji,2) - zcoefs
267            ENDIF
268         END DO
269      ENDIF
270!!cr      ENDIF
271#endif
272     
273      ! North fold boundary condition
274      !     all the coef are already set to zero as bmask is initialized to
275      !     zero on duplicated lignes and portion of lignes
276     
277      ! 3. Preconditioned matrix
278      ! ------------------------
279     
280      IF( nsolv /= 3 ) THEN
281         
282         ! SOR and PCG solvers
283         DO jj = 1, jpj
284            DO ji = 1, jpi
285               IF( bmask(ji,jj) /= 0.e0 )   gcdprc(ji,jj) = 1.e0 / gcdmat(ji,jj)
286            END DO
287         END DO
288         
289         gcp(:,:,1) = gcp(:,:,1) * gcdprc(:,:)
290         gcp(:,:,2) = gcp(:,:,2) * gcdprc(:,:)
291         gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:)
292         gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:)
293         IF( nsolv == 2 )  gccd(:,:) = sor * gcp(:,:,2)
294
295         IF( nsolv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0) THEN
296            CALL lbc_lnk_e( gcp   (:,:,1), c_solver_pt, 1. )   ! lateral boundary conditions
297            CALL lbc_lnk_e( gcp   (:,:,2), c_solver_pt, 1. )   ! lateral boundary conditions
298            CALL lbc_lnk_e( gcp   (:,:,3), c_solver_pt, 1. )   ! lateral boundary conditions
299            CALL lbc_lnk_e( gcp   (:,:,4), c_solver_pt, 1. )   ! lateral boundary conditions
300            CALL lbc_lnk_e( gcdprc(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions
301            CALL lbc_lnk_e( gcdmat(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions         
302            IF( npolj /= 0 ) CALL sol_exd( gcp , c_solver_pt ) ! switch northernelements
303         END IF
304
305      ELSE
306         
307         ! FETI method
308         ! if feti solver : gcdprc is a mask for the non-overlapping
309         !   data structuring
310         
311         DO jj = 1, jpj
312            DO ji = 1, jpi
313               IF( bmask(ji,jj) /= 0.e0 ) THEN
314                  gcdprc(ji,jj) = 1.e0
315               ELSE
316                  gcdprc(ji,jj) = 0.e0
317               ENDIF
318            END DO
319         END DO
320         
321         ! so "common" line & "common" column have to be !=0 except on global
322         !   domain boundaries
323         ! pbs with nbondi if nperio != 2 ?
324         !   ii = nldi-1
325         ! pb with nldi value if jperio==1 : nbondi modifyed at the end
326         !   of inimpp.F => pb
327         ! pb with periodicity conditions : iiend, ijend
328         
329         ijend = nlej
330         iiend = nlei
331         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 )   iiend = nlci - jpreci
332         ii = jpreci
333         
334         ! case number 1
335         
336         IF( nbondi /= -1 .AND. nbondi /= 2 ) THEN
337            DO jj = 1, ijend
338               IF( fmask(ii,jj,1) == 1. ) gcdprc(ii,jj) = 1.
339            END DO
340         ENDIF
341         
342         ! case number 2
343         
344         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
345            DO jj = 1, ijend
346               IF( fmask(ii,jj,1) == 1. ) gcdprc(ii,jj) = 1.
347            END DO
348         ENDIF
349         
350         !      ij = nldj-1
351         ! pb with nldi value if jperio==1 : nbondi modifyed at the end
352         !   of inimpp.F => pb, here homogeneisation...
353         
354         ij = jprecj
355         IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN
356            DO ji = 1, iiend
357               IF( fmask(ji,ij,1) == 1. ) gcdprc(ji,ij) = 1.
358            END DO
359         ENDIF
360      ENDIF
361     
362     
363      ! 4. Initialization the arrays used in pcg
364      ! ----------------------------------------
365      gcb  (:,:) = 0.e0
366      gcr  (:,:) = 0.e0
367      gcdes(:,:) = 0.e0
368      gccd (:,:) = 0.e0
369     
370      ! FETI method
371      IF( nsolv == 3 ) THEN
372         CALL fetmat       ! Matrix treatment : Neumann condition, inverse computation
373         CALL fetsch       ! data framework for the Schur Dual solver
374      ENDIF
375     
376   END SUBROUTINE sol_mat
377
378
379   SUBROUTINE sol_exd( pt3d, cd_type )
380      !!----------------------------------------------------------------------
381      !!                  ***  routine sol_exd  ***
382      !!                 
383      !! ** Purpose :   Reorder gcb coefficient on the extra outer  halo
384      !!                at north fold in case of T or F pivot
385      !!
386      !! ** Method  :   Perform a circular permutation of the coefficients on
387      !!                the total area strictly above the pivot point,
388      !!                and on the semi-row of the pivot point   
389      !!               
390      !! History :
391      !!   9.0  !  05-09  (R. Benshila)  original routine
392      !!----------------------------------------------------------------------
393      !! * Arguments
394      CHARACTER(len=1) , INTENT( in ) ::   &
395         cd_type       ! define the nature of pt2d array grid-points
396         !             !  = T , U , V , F , W
397         !             !  = S : T-point, north fold treatment
398         !             !  = G : F-point, north fold treatment
399         !             !  = I : sea-ice velocity at F-point with index shift
400      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4), INTENT( inout ) ::   &
401         pt3d          ! 2D array on which the boundary condition is applied
402
403      !! * Local variables
404      INTEGER  ::   ji, jk      ! dummy loop indices
405      INTEGER  ::   iloc                ! temporary integers
406      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4) ::   &
407         ztab          ! 2D array on which the boundary condition is applied
408      !!----------------------------------------------------------------------
409
410      ztab = pt3d
411
412      ! north fold treatment
413      ! -----------------------
414 
415      SELECT CASE ( npolj )
416         
417         CASE ( 3 , 4 )   !  T pivot
418         iloc = jpiglo/2 +1 
419           
420            SELECT CASE ( cd_type )
421 
422            CASE ( 'T', 'S', 'U', 'W' )
423               DO jk =1, 4
424                  DO ji = 1-jpr2di, nlci+jpr2di
425                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
426                  ENDDO
427               ENDDO
428
429              DO jk =1, 4
430                  DO ji = nlci+jpr2di, 1-jpr2di,  -1
431                     IF( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) &
432                       & .OR. ( mi0(iloc) == jpi+1 ) ) EXIT
433                     pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
434                  ENDDO
435               ENDDO
436
437            CASE ( 'F' ,'G' , 'I', 'V' )
438               DO jk =1, 4
439                  DO ji = 1-jpr2di, nlci+jpr2di
440                     pt3d(ji,nlcj-1:nlcj+jpr2dj,jk) = ztab(ji,nlcj-1:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
441                  ENDDO
442               ENDDO
443
444            END SELECT   ! cd_type
445 
446         CASE ( 5 , 6 )                 ! F pivot
447          iloc=jpiglo/2
448
449            SELECT CASE (cd_type )
450
451            CASE ( 'T'  ,'S', 'U', 'W')
452               DO jk =1, 4
453                  DO ji = 1-jpr2di, nlci+jpr2di
454                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
455                  ENDDO
456               ENDDO
457
458            CASE ( 'F' ,'G' , 'I', 'V' )
459               DO jk =1, 4
460                  DO ji = 1-jpr2di, nlci+jpr2di
461                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
462                  ENDDO
463               ENDDO
464               DO jk =1, 4
465                  DO ji = nlci+jpr2di, 1-jpr2di,  -1
466                    IF ( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) &
467                       & .OR. ( mi0(iloc) == jpi+1 ) ) EXIT
468                    pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
469                  ENDDO
470               ENDDO
471
472            END SELECT   ! cd_type
473
474         END SELECT   ! npolj
475 
476   END SUBROUTINE sol_exd
477
478#if defined key_feti
479
480   SUBROUTINE fetstr
481      !!---------------------------------------------------------------------
482      !!                  ***  ROUTINE fetstr  ***
483      !!               
484      !! ** Purpose :   Construction of the matrix of the barotropic stream
485      !!      function system.
486      !!      Finite Elements Tearing & Interconnecting (FETI) approach
487      !!      Memory allocation and interface definition for the solver
488      !!
489      !! ** Method :
490      !!
491      !! References :
492      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
493      !!      A domain decomposition solver to compute the barotropic
494      !!      component of an OGCM in the parallel processing field.
495      !!      Ocean Modelling, issue 105, december 94.
496      !!
497      !! History :
498      !!        !  98-02 (M. Guyon)  Original code
499      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
500      !!----------------------------------------------------------------------
501      !! * Modules used
502      USE lib_feti            ! feti librairy
503      !! * Local declarations
504      INTEGER ::   iiend, ijend, iperio   ! temporary integers
505      !!---------------------------------------------------------------------
506     
507     
508      ! Preconditioning technics of the Dual Schur Operator
509      ! <= definition of the Coarse Grid solver
510      ! <= dimension of the nullspace of the local operators
511      ! <= Neumann boundaries conditions
512     
513      ! 0. Initializations
514      ! ------------------
515     
516      ndkerep = 1
517
518      ! initialization of the superstructures management
519
520      malxm = 1
521      malim = 1
522
523      ! memory space for the pcpg associated with the FETI dual formulation
524      ! ndkerep is associated to the list of rigid modes,
525      ! ndkerep == 1 because the Dual Operator
526      ! is a first order operator due to SPG elliptic Operator is a
527      ! second order operator
528
529      nim = 50
530      nim = nim + ndkerep
531      nim = nim + 2*jpi + 2*jpj
532      nim = nim + jpi*jpj
533
534      nxm = 33
535      nxm = nxm + 4*jpnij
536      nxm = nxm + 19*(jpi+jpj)
537      nxm = nxm + 13*jpi*jpj
538      nxm = nxm + jpi*jpi*jpj
539
540      ! krylov space memory
541
542      iperio = 0
543      IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6) iperio = 1
544      nxm = nxm + 3*(jpnij-jpni)*jpi
545      nxm = nxm + 3*(jpnij-jpnj+iperio)*jpj
546      nxm = nxm + 2*(jpi+jpj)*(jpnij-jpni)*jpi
547      nxm = nxm + 2*(jpi+jpj)*(jpnij-jpnj+iperio)*jpj
548
549      ! Resolution with the Schur dual Method ( frontal and local solver by
550      ! blocks
551      ! Case with a local symetrical matrix
552      ! The local matrix is stored in a multi-column form
553      ! The total number of nodes for this subdomain is named "noeuds"
554
555      noeuds = jpi*jpj
556      nifmat = jpi-1
557      njfmat = jpj-1
558      nelem = nifmat*njfmat
559      npe = 4
560      nmorse = 5*noeuds
561     
562      ! 1. mesh building
563      ! ----------------
564     
565      ! definition of specific information for a subdomain
566      !  narea           : subdomain number = processor number +1
567      !  ninterf         : neighbour subdomain number
568      !  nni             : interface point number
569      !  ndvois array    : neighbour subdomain list
570      !  maplistin array : node pointer at each interface
571      !  maplistin array : concatened list of interface nodes
572     
573      !  messag coding is necessary by interface type for avoid collision
574      !  if nperio == 1
575     
576      !  lint  array     : indoor interface list / type
577      !  lext  array     : outdoor interface list / type
578     
579      !  domain with jpniXjpnj subdomains
580     
581      CALL feti_inisub(nifmat,njfmat,nbondi,nbondj,nperio,   &
582          nbsw,nbnw,nbse,nbne,ninterf,ninterfc,nni,nnic)
583
584      CALL feti_creadr(malim,malimax,nim,3*ninterf ,mandvois ,'ndvois' )
585      CALL feti_creadr(malim,malimax,nim,3*ninterfc,mandvoisc,'ndvoisc')
586      CALL feti_creadr(malim,malimax,nim,ninterfc+1,maplistin,'plistin')
587      CALL feti_creadr(malim,malimax,nim,nnic      ,malistin ,'listin' )
588
589      ! pb with periodicity conditions : iiend, ijend
590
591      ijend = nlej
592      iiend = nlei
593      IF (jperio == 1) iiend = nlci - jpreci
594
595      CALL feti_subound(nifmat,njfmat,nldi,iiend,nldj,ijend,   &
596          narea,nbondi,nbondj,nperio,   &
597          ninterf,ninterfc,   &
598          nowe,noea,noso,nono,   &
599          nbsw,nbnw,nbse,nbne,   &
600          npsw,npnw,npse,npne,   &
601          mfet(mandvois),mfet(mandvoisc),   &
602          mfet(maplistin),nnic,mfet(malistin) )
603
604   END SUBROUTINE fetstr
605
606
607   SUBROUTINE fetmat
608      !!---------------------------------------------------------------------
609      !!                   ***  ROUTINE fetmat  ***
610      !!           
611      !! ** Purpose :   Construction of the matrix of the barotropic stream
612      !!      function system.
613      !!      Finite Elements Tearing & Interconnecting (FETI) approach
614      !!      Matrix treatment : Neumann condition, inverse computation
615      !!
616      !! ** Method :
617      !!
618      !! References :
619      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
620      !!      A domain decomposition solver to compute the barotropic
621      !!      component of an OGCM in the parallel processing field.
622      !!      Ocean Modelling, issue 105, december 94.
623      !!
624      !! History :
625      !!        !  98-02 (M. Guyon)  Original code
626      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
627      !!----------------------------------------------------------------------
628      !! * Modules used
629      USE lib_feti            ! feti librairy
630      !! * Local declarations
631      INTEGER ::   ji, jj, jk, jl
632      INTEGER ::   iimask(jpi,jpj)
633      INTEGER ::   iiend, ijend
634      REAL(wp) ::   zres, zres2, zdemi
635      !!---------------------------------------------------------------------
636
637      ! Matrix computation
638      ! ------------------
639
640      CALL feti_creadr(malxm,malxmax,nxm,nmorse,maan,'matrice a')
641
642      nnitot = nni
643
644      CALL mpp_sum( nnitot, 1, numit0ete )
645      CALL feti_creadr(malxm,malxmax,nxm,npe*npe,maae,'ae')
646
647      ! initialisation of the local barotropic matrix
648      ! local boundary conditions on the halo
649
650      CALL lbc_lnk( gcp(:,:,1), 'F', 1)
651      CALL lbc_lnk( gcp(:,:,2), 'F', 1) 
652      CALL lbc_lnk( gcp(:,:,3), 'F', 1) 
653      CALL lbc_lnk( gcp(:,:,4), 'F', 1) 
654      CALL lbc_lnk( gcdmat    , 'T', 1)
655
656      ! Neumann conditions
657      ! initialisation of the integer Neumann Mask
658
659      CALL feti_iclr(jpi*jpj,iimask)
660      DO jj = 1, jpj
661         DO ji = 1, jpi
662            iimask(ji,jj) = INT( gcdprc(ji,jj) )
663         END DO
664      END DO
665
666      ! regularization of the local matrix
667
668      DO jj = 1, jpj
669         DO ji = 1, jpi
670            gcdmat(ji,jj) = gcdmat(ji,jj) * gcdprc(ji,jj) + 1. - gcdprc(ji,jj)
671         END DO
672      END DO
673
674      DO jk = 1, 4
675         DO jj = 1, jpj
676            DO ji = 1, jpi
677               gcp(ji,jj,jk) = gcp(ji,jj,jk) * gcdprc(ji,jj)
678            END DO
679         END DO
680      END DO
681     
682      ! implementation of the west, east, north & south Neumann conditions
683
684      zdemi  = 0.5
685
686      ! pb with periodicity conditions : iiend, ijend
687
688      ijend = nlej
689      iiend = nlei
690      IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) iiend = nlci - jpreci
691
692      IF( nbondi == 2 .AND. (nperio /= 1 .OR. nperio /= 4 .OR. nperio == 6) ) THEN
693
694         ! with the periodicity : no east/west interface if nbondi = 2
695         ! and nperio != 1
696
697      ELSE 
698         ! west
699         IF( nbondi /= -1 ) THEN
700            DO jj = 1, jpj
701               IF( iimask(1,jj) /= 0 ) THEN
702                  gcp(1,jj,2) = 0.e0
703                  gcp(1,jj,1) = zdemi * gcp(1,jj,1)
704                  gcp(1,jj,4) = zdemi * gcp(1,jj,4)
705               ENDIF
706            END DO
707            DO jj = 1, jpj
708               IF( iimask(1,jj) /= 0 ) THEN
709                  gcdmat(1,jj) = - ( gcp(1,jj,1) + gcp(1,jj,2) + gcp(1,jj,3) + gcp(1,jj,4) )
710               ENDIF
711            END DO
712         ENDIF
713         ! east
714         IF( nbondi /= 1 ) THEN
715            DO jj = 1, jpj
716               IF( iimask(iiend,jj) /= 0 ) THEN
717                  gcp(iiend,jj,3) = 0.e0
718                  gcp(iiend,jj,1) = zdemi * gcp(iiend,jj,1)
719                  gcp(iiend,jj,4) = zdemi * gcp(iiend,jj,4)
720               ENDIF
721            END DO
722            DO jj = 1, jpj
723               IF( iimask(iiend,jj) /= 0 ) THEN
724                  gcdmat(iiend,jj) = - ( gcp(iiend,jj,1) + gcp(iiend,jj,2)   &
725                                       + gcp(iiend,jj,3) + gcp(iiend,jj,4) )
726               ENDIF
727            END DO
728         ENDIF
729      ENDIF
730
731      ! south
732      IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN
733         DO ji = 1, jpi
734            IF( iimask(ji,1) /= 0 ) THEN
735               gcp(ji,1,1) = 0.e0
736               gcp(ji,1,2) = zdemi * gcp(ji,1,2)
737               gcp(ji,1,3) = zdemi * gcp(ji,1,3)
738            ENDIF
739         END DO
740         DO ji = 1, jpi
741            IF( iimask(ji,1) /= 0 ) THEN
742               gcdmat(ji,1) = - ( gcp(ji,1,1) + gcp(ji,1,2) + gcp(ji,1,3) + gcp(ji,1,4) )
743            ENDIF
744         END DO
745      ENDIF
746     
747      ! north
748      IF( nbondj /= 1 .AND. nbondj /= 2 ) THEN
749         DO ji = 1, jpi
750            IF( iimask(ji,ijend) /= 0 ) THEN
751               gcp(ji,ijend,4) = 0.e0
752               gcp(ji,ijend,2) = zdemi * gcp(ji,ijend,2) 
753               gcp(ji,ijend,3) = zdemi * gcp(ji,ijend,3)
754            ENDIF
755         END DO
756         DO ji = 1, jpi
757            IF( iimask(ji,ijend) /= 0 ) THEN
758               gcdmat(ji,ijend) = - ( gcp(ji,ijend,1) + gcp(ji,ijend,2)   &
759                                    + gcp(ji,ijend,3) + gcp(ji,ijend,4) )
760            ENDIF
761         END DO
762      ENDIF
763
764      ! matrix terms are  saved in FETI solver arrays
765      CALL feti_vmov(noeuds,gcp(1,1,1),wfeti(maan))
766      CALL feti_vmov(noeuds,gcp(1,1,2),wfeti(maan+noeuds))
767      CALL feti_vmov(noeuds,gcdmat,wfeti(maan+2*noeuds))
768      CALL feti_vmov(noeuds,gcp(1,1,3),wfeti(maan+3*noeuds))
769      CALL feti_vmov(noeuds,gcp(1,1,4),wfeti(maan+4*noeuds))
770
771      ! construction of Dirichlet liberty degrees array
772      CALL feti_subdir(nifmat,njfmat,noeuds,ndir,iimask)
773      CALL feti_creadr(malim,malimax,nim,ndir,malisdir,'lisdir')
774      CALL feti_listdir(jpi,jpj,iimask,ndir,mfet(malisdir))
775
776      ! stop onto  matrix term for Dirichlet conditions
777      CALL feti_blomat(nifmat+1,njfmat+1,wfeti(maan),ndir,mfet(malisdir))
778
779      ! reservation of factorized diagonal blocs and temporary array for
780      ! factorization
781      npblo = (njfmat+1) * (nifmat+1) * (nifmat+1)
782      ndimax = nifmat+1
783
784      CALL feti_creadr(malxm,malxmax,nxm,npblo,mablo,'blo')
785      CALL feti_creadr(malxm,malxmax,nxm,noeuds,madia,'dia')
786      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mav,'v')
787      CALL feti_creadr(malxm,malxmax,nxm,ndimax*ndimax,mautil,'util')
788
789      ! stoping the rigid modes
790
791      ! the number of rigid modes =< Max [dim(Ker(Ep))]
792      !                                p=1,Np
793
794      CALL feti_creadr(malim,malimax,nim,ndkerep,malisblo,'lisblo')
795
796      ! Matrix factorization
797
798      CALL feti_front(noeuds,nifmat+1,njfmat+1,wfeti(maan),npblo,   &
799          wfeti(mablo),wfeti(madia),   &
800          wfeti(mautil),wfeti(mav),ndlblo,mfet(malisblo),ndkerep)
801      CALL feti_prext(noeuds,wfeti(madia))
802
803      ! virtual dealloc => we have to see for a light f90 version
804      ! the super structure is removed to clean the coarse grid
805      ! solver structure
806
807      malxm = madia
808      CALL feti_vclr(noeuds,wfeti(madia))
809      CALL feti_vclr(noeuds,wfeti(mav))
810      CALL feti_vclr(ndimax*ndimax,wfeti(mautil))
811
812      ! ndlblo is the dimension of the local nullspace .=<. the size of the
813      ! memory of the superstructure associated to the nullspace : ndkerep
814      ! ndkerep is introduced to avoid messages "out of bounds" when memory
815      ! is checked
816
817      ! copy matrix for Dirichlet condition
818
819      CALL feti_creadr(malxm,malxmax,nxm,noeuds,miax,'x')
820      CALL feti_creadr(malxm,malxmax,nxm,noeuds,may,'y')
821      CALL feti_creadr(malxm,malxmax,nxm,noeuds,maz,'z')
822
823      ! stoping the rigid modes
824
825      ! ndlblo is the dimension of the local nullspace .=<. the size of the
826      ! memory of the superstructure associated to the nullspace : ndkerep
827      ! ndkerep is introduced to avoid messages "out of bounds" when memory
828      ! is checked
829
830      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*noeuds,mansp,'nsp')
831      CALL feti_blomat1(nifmat+1,njfmat+1,wfeti(maan),ndlblo,   &
832          mfet(malisblo),wfeti(mansp))     
833
834      ! computation of operator kernel
835
836      CALL feti_nullsp(noeuds,nifmat+1,njfmat+1,npblo,wfeti(mablo),   &
837          wfeti(maan),ndlblo,mfet(malisblo),wfeti(mansp),   &
838          wfeti(maz))
839
840      ! test of the factorisation onto each sub domain
841
842      CALL feti_init(noeuds,wfeti(may))
843      CALL feti_blodir(noeuds,wfeti(may),ndir,mfet(malisdir))
844      CALL feti_blodir(noeuds,wfeti(may),ndlblo,mfet(malisblo))
845      CALL feti_vclr(noeuds,wfeti(miax))
846      CALL feti_resloc(noeuds,nifmat+1,njfmat+1,wfeti(maan),npblo,   &
847          wfeti(mablo),wfeti(may),wfeti(miax),wfeti(maz)) 
848      CALL feti_proax(noeuds,nifmat+1,njfmat+1,wfeti(maan),wfeti(miax),   &
849          wfeti(maz))
850      CALL feti_blodir(noeuds,wfeti(maz),ndlblo,mfet(malisblo))
851      CALL feti_vsub(noeuds,wfeti(may),wfeti(maz),wfeti(maz))
852
853      zres2 = 0.e0
854      DO jl = 1, noeuds
855         zres2 = zres2 + wfeti(may+jl-1) * wfeti(may+jl-1)
856      END DO
857      CALL mpp_sum(zres2,1,zres)
858
859      res2 = 0.e0
860      DO jl = 1, noeuds
861         res2 = res2 + wfeti(maz+jl-1) * wfeti(maz+jl-1)
862      END DO
863      res2 = res2 / zres2
864      CALL mpp_sum(res2,1,zres)
865
866      res2 = SQRT(res2)
867      IF(lwp) WRITE(numout,*) 'global residu : sqrt((Ax-b,Ax-b)/(b.b)) =', res2
868
869      IF( res2 > (eps/100.) ) THEN
870         IF(lwp) WRITE (numout,*) 'eps is :',eps
871         IF(lwp) WRITE (numout,*) 'factorized matrix precision :',res2
872         STOP
873      ENDIF
874
875   END SUBROUTINE fetmat
876
877
878   SUBROUTINE fetsch
879      !!---------------------------------------------------------------------
880      !!                  ***  ROUTINE fetsch  ***
881      !!       
882      !! ** Purpose :
883      !!     Construction of the matrix of the barotropic stream function
884      !!     system.
885      !!     Finite Elements Tearing & Interconnecting (FETI) approach
886      !!     Data framework for the Schur Dual solve
887      !!
888      !! ** Method :
889      !!
890      !! References :
891      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
892      !!      A domain decomposition solver to compute the barotropic
893      !!      component of an OGCM in the parallel processing field.
894      !!      Ocean Modelling, issue 105, december 94.
895      !!
896      !! History :
897      !!        !  98-02 (M. Guyon)  Original code
898      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
899      !!----------------------------------------------------------------------
900      !! * Modules used
901      USE lib_feti            ! feti librairy
902      !! * Local declarations
903      !!---------------------------------------------------------------------
904
905      !  computing weights for the conform construction
906
907      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mapoids ,'poids' )
908      CALL feti_creadr(malxm,malxmax,nxm,nnic  ,mabufin ,'bufin' )
909      CALL feti_creadr(malxm,malxmax,nxm,nnic  ,mabufout,'bufout')
910
911!!    CALL feti_poids(ninterfc,mfet(mandvoisc),mfet(maplistin),nnic,   &
912!!        mfet(malistin),narea,noeuds,wfeti(mapoids),wfeti(mabufin),   &
913!!        wfeti(mabufout) )
914      CALL feti_poids(ninterfc,                                nnic,   &
915          mfet(malistin),      noeuds,wfeti(mapoids)                )
916
917
918      ! Schur dual arrays
919 
920      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mabitw,'bitw') 
921      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mautilu,'utilu') 
922      CALL feti_creadr(malxm,malxmax,nxm,nni,malambda,'lambda') 
923      CALL feti_creadr(malxm,malxmax,nxm,nni,mag,'g') 
924      CALL feti_creadr(malxm,malxmax,nxm,nni,mapg,'pg') 
925      CALL feti_creadr(malxm,malxmax,nxm,nni,mamg,'mg') 
926      CALL feti_creadr(malxm,malxmax,nxm,nni,maw,'w') 
927      CALL feti_creadr(malxm,malxmax,nxm,nni,madw,'dw')
928
929      !  coarse grid solver dimension and arrays
930
931      nitmaxete = ndlblo
932      CALL  mpp_sum(nitmaxete,1,numit0ete)
933
934      nitmaxete = nitmaxete + 1
935      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maxnul,'xnul')
936      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maynul,'ynul')
937      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maeteg,'eteg')
938      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maeteag,'eteag')
939      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*nitmaxete,maeted,'eted')
940      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*nitmaxete,maetead,'etead')
941      CALL feti_creadr(malxm,malxmax,nxm,nitmaxete,maeteadd,'eteadd')
942      CALL feti_creadr(malxm,malxmax,nxm,nitmaxete,maetegamm,'etegamm')
943      CALL feti_creadr(malxm,malxmax,nxm,nni,maetew,'etew') 
944      CALL feti_creadr(malxm,malxmax,nxm,noeuds,maetev,'etev') 
945
946      ! construction of semi interface arrays
947
948      CALL feti_creadr(malim,malimax,nim,ninterf+1,maplistih,'plistih')
949!!    CALL feti_halfint(ninterf,mfet(mandvois),mfet(maplistin),nni,   &
950!!        mfet(maplistih),nnih,narea)
951      CALL feti_halfint(ninterf,mfet(mandvois),mfet(maplistin),       & 
952          mfet(maplistih),nnih      )
953
954      CALL feti_creadr(malxm,malxmax,nxm,nnih,magh,'gh') 
955
956      ! Schur Dual Method
957
958      nmaxd = nnitot / 2
959
960      !  computation of the remain array for descent directions
961
962      nmaxd = min0(nmaxd,(nxm-nitmaxete-malxm)/(2*nnih+3))
963      CALL mpp_min(nmaxd,1,numit0ete)
964
965      nitmax = nnitot/2
966      epsilo = eps
967      ntest = 0
968
969      ! Krylov space construction
970
971      CALL feti_creadr(malxm,malxmax,nxm,nnih*nmaxd,mawj,'wj') 
972      CALL feti_creadr(malxm,malxmax,nxm,nnih*nmaxd,madwj,'dwj') 
973      CALL feti_creadr(malxm,malxmax,nxm,nmaxd,madwwj,'dwwj') 
974      CALL feti_creadr(malxm,malxmax,nxm,nmaxd,magamm,'gamm') 
975      CALL feti_creadr(malxm,malxmax,nxm,max0(nmaxd,nitmaxete),mawork,'work')
976      mjj0 = 0 
977      numit0ete = 0 
978
979   END SUBROUTINE fetsch
980
981#else
982   SUBROUTINE fetstr                 ! Empty routine
983   END SUBROUTINE fetstr
984   SUBROUTINE fetmat                 ! Empty routine
985   END SUBROUTINE fetmat
986   SUBROUTINE fetsch                 ! Empty routine
987   END SUBROUTINE fetsch
988#endif
989
990   !!======================================================================
991END MODULE solmat
Note: See TracBrowser for help on using the repository browser.