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

Last change on this file since 392 was 392, checked in by opalod, 15 years ago

RB:nemo_v1_update_038: first integration of Agrif :

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