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 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 29.5 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 lib_mpp         ! distributed memory computing
20   USE dynspg_rl
21   USE dynspg_fsc
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC sol_mat     ! routine called by inisol.F90
28   !!----------------------------------------------------------------------
29   !!   OPA 9.0 , LOCEAN-IPSL (2005)
30   !! $Header$
31   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
32   !!----------------------------------------------------------------------
33
34CONTAINS
35
36   SUBROUTINE sol_mat
37      !!----------------------------------------------------------------------
38      !!                  ***  ROUTINE sol_mat  ***
39      !!
40      !! ** Purpose :   Construction of the matrix of used by the elliptic
41      !!      solvers (either sor, pcg or feti methods).
42      !!
43      !! ** Method  :   The matrix depends on the type of free surface:
44      !!       * lk_dynspg_rl=T: rigid lid formulation
45      !!      The matrix is built for the barotropic stream function system.
46      !!      a diagonal preconditioning matrix is also defined.
47      !!       * lk_dynspg_fsc=T: free surface formulation
48      !!      The matrix is built for the divergence of the transport system
49      !!      a diagonal preconditioning matrix is also defined.
50      !!        Note that for feti solver (nsolv=3) a specific initialization
51      !!      is required (call to fetstr.F) for memory allocation and inter-
52      !!      face definition.
53      !!
54      !! ** Action  : - gcp    : extra-diagonal elements of the matrix
55      !!              - gcdmat : preconditioning matrix (diagonal elements)
56      !!              - gcdprc : inverse of the preconditioning matrix
57      !!
58      !! History :
59      !!   1.0  !  88-04  (G. Madec)  Original code
60      !!        !  91-11  (G. Madec)
61      !!        !  93-03  (M. Guyon)  symetrical conditions
62      !!        !  93-06  (M. Guyon)  suppress pointers
63      !!        !  96-05  (G. Madec)  merge sor and pcg formulations
64      !!        !  96-11  (A. Weaver)  correction to preconditioning
65      !!        !  98-02  (M. Guyon)  FETI method
66      !!   8.5  !  02-08  (G. Madec)  F90: Free form
67      !!        !  02-11  (C. Talandier, A-M. Treguier) Free surface & Open boundaries
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_fsc && ! defined key_obc
98!!cr      IF( lk_dynspg_fsc .AND. .NOT.lk_obc ) THEN   !bug missing lk_dynspg_fsc_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_fsc && defined key_obc
119!!cr      ELSEIF( lk_dynspg_fsc .AND. lk_obc ) THEN     !bug missing lk_dynspg_fsc_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
186      ! 2. Boundary conditions
187      ! ----------------------
188     
189      ! Cyclic east-west boundary conditions
190      !     ji=2 is the column east of ji=jpim1 and reciprocally,
191      !     ji=jpim1 is the column west of ji=2
192      !     all the coef are already set to zero as bmask is initialized to
193      !     zero for ji=1 and ji=jpj in dommsk.
194     
195      ! Symetrical conditions
196      ! free surface: no specific action
197      ! bsf system: n-s gradient of bsf = 0 along j=2 (perhaps a bug !!!!!!)
198      ! the diagonal coefficient of the southern grid points must be modify to
199      ! account for the existence of the south symmetric bassin.
200     
201!!cr      IF( .NOT.lk_dynspg_fsc ) THEN   !bug missing lk_dynspg_fsc_atsk
202#if ! defined key_dynspg_fsc
203      IF( nperio == 2 ) THEN
204         DO ji = 1, jpi
205            IF( bmask(ji,2) /= 0.e0 ) THEN
206               zcoefs = - hur(ji,2)*e1u(ji,2)/e2u(ji,2)
207               gcdmat(ji,2) = gcdmat(ji,2) - zcoefs
208            ENDIF
209         END DO
210      ENDIF
211!!cr      ENDIF
212#endif
213     
214      ! North fold boundary condition
215      !     all the coef are already set to zero as bmask is initialized to
216      !     zero on duplicated lignes and portion of lignes
217     
218      ! 3. Preconditioned matrix
219      ! ------------------------
220     
221      IF( nsolv /= 3 ) THEN
222         
223         ! SOR and PCG solvers
224         DO jj = 1, jpj
225            DO ji = 1, jpi
226               IF( bmask(ji,jj) /= 0.e0 )   gcdprc(ji,jj) = 1.e0 / gcdmat(ji,jj)
227            END DO
228         END DO
229         
230         gcp(:,:,1) = gcp(:,:,1) * gcdprc(:,:)
231         gcp(:,:,2) = gcp(:,:,2) * gcdprc(:,:)
232         gcp(:,:,3) = gcp(:,:,3) * gcdprc(:,:)
233         gcp(:,:,4) = gcp(:,:,4) * gcdprc(:,:)
234         IF( nsolv == 2 )   gccd(:,:) = sor * gcp(:,:,2)
235
236      ELSE
237         
238         ! FETI method
239         ! if feti solver : gcdprc is a mask for the non-overlapping
240         !   data structuring
241         
242         DO jj = 1, jpj
243            DO ji = 1, jpi
244               IF( bmask(ji,jj) /= 0.e0 ) THEN
245                  gcdprc(ji,jj) = 1.e0
246               ELSE
247                  gcdprc(ji,jj) = 0.e0
248               ENDIF
249            END DO
250         END DO
251         
252         ! so "common" line & "common" column have to be !=0 except on global
253         !   domain boundaries
254         ! pbs with nbondi if nperio != 2 ?
255         !   ii = nldi-1
256         ! pb with nldi value if jperio==1 : nbondi modifyed at the end
257         !   of inimpp.F => pb
258         ! pb with periodicity conditions : iiend, ijend
259         
260         ijend = nlej
261         iiend = nlei
262         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 )   iiend = nlci - jpreci
263         ii = jpreci
264         
265         ! case number 1
266         
267         IF( nbondi /= -1 .AND. nbondi /= 2 ) THEN
268            DO jj = 1, ijend
269               IF( fmask(ii,jj,1) == 1. ) gcdprc(ii,jj) = 1.
270            END DO
271         ENDIF
272         
273         ! case number 2
274         
275         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
276            DO jj = 1, ijend
277               IF( fmask(ii,jj,1) == 1. ) gcdprc(ii,jj) = 1.
278            END DO
279         ENDIF
280         
281         !      ij = nldj-1
282         ! pb with nldi value if jperio==1 : nbondi modifyed at the end
283         !   of inimpp.F => pb, here homogeneisation...
284         
285         ij = jprecj
286         IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN
287            DO ji = 1, iiend
288               IF( fmask(ji,ij,1) == 1. ) gcdprc(ji,ij) = 1.
289            END DO
290         ENDIF
291      ENDIF
292     
293     
294      ! 4. Initialization the arrays used in pcg
295      ! ----------------------------------------
296      gcx  (:,:) = 0.e0
297      gcxb (:,:) = 0.e0
298      gcb  (:,:) = 0.e0
299      gcr  (:,:) = 0.e0
300      gcdes(:,:) = 0.e0
301      gccd (:,:) = 0.e0
302     
303      ! FETI method
304      IF( nsolv == 3 ) THEN
305         CALL fetmat       ! Matrix treatment : Neumann condition, inverse computation
306         CALL fetsch       ! data framework for the Schur Dual solver
307      ENDIF
308     
309   END SUBROUTINE sol_mat
310
311#if defined key_feti
312
313   SUBROUTINE fetstr
314      !!---------------------------------------------------------------------
315      !!                  ***  ROUTINE fetstr  ***
316      !!               
317      !! ** Purpose :   Construction of the matrix of the barotropic stream
318      !!      function system.
319      !!      Finite Elements Tearing & Interconnecting (FETI) approach
320      !!      Memory allocation and interface definition for the solver
321      !!
322      !! ** Method :
323      !!
324      !! References :
325      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
326      !!      A domain decomposition solver to compute the barotropic
327      !!      component of an OGCM in the parallel processing field.
328      !!      Ocean Modelling, issue 105, december 94.
329      !!
330      !! History :
331      !!        !  98-02 (M. Guyon)  Original code
332      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
333      !!----------------------------------------------------------------------
334      !! * Modules used
335      USE lib_feti            ! feti librairy
336      !! * Local declarations
337      INTEGER ::   iiend, ijend, iperio   ! temporary integers
338      !!---------------------------------------------------------------------
339     
340     
341      ! Preconditioning technics of the Dual Schur Operator
342      ! <= definition of the Coarse Grid solver
343      ! <= dimension of the nullspace of the local operators
344      ! <= Neumann boundaries conditions
345     
346      ! 0. Initializations
347      ! ------------------
348     
349      ndkerep = 1
350
351      ! initialization of the superstructures management
352
353      malxm = 1
354      malim = 1
355
356      ! memory space for the pcpg associated with the FETI dual formulation
357      ! ndkerep is associated to the list of rigid modes,
358      ! ndkerep == 1 because the Dual Operator
359      ! is a first order operator due to SPG elliptic Operator is a
360      ! second order operator
361
362      nim = 50
363      nim = nim + ndkerep
364      nim = nim + 2*jpi + 2*jpj
365      nim = nim + jpi*jpj
366
367      nxm = 33
368      nxm = nxm + 4*jpnij
369      nxm = nxm + 19*(jpi+jpj)
370      nxm = nxm + 13*jpi*jpj
371      nxm = nxm + jpi*jpi*jpj
372
373      ! krylov space memory
374
375      iperio = 0
376      IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6) iperio = 1
377      nxm = nxm + 3*(jpnij-jpni)*jpi
378      nxm = nxm + 3*(jpnij-jpnj+iperio)*jpj
379      nxm = nxm + 2*(jpi+jpj)*(jpnij-jpni)*jpi
380      nxm = nxm + 2*(jpi+jpj)*(jpnij-jpnj+iperio)*jpj
381
382      ! Resolution with the Schur dual Method ( frontal and local solver by
383      ! blocks
384      ! Case with a local symetrical matrix
385      ! The local matrix is stored in a multi-column form
386      ! The total number of nodes for this subdomain is named "noeuds"
387
388      noeuds = jpi*jpj
389      nifmat = jpi-1
390      njfmat = jpj-1
391      nelem = nifmat*njfmat
392      npe = 4
393      nmorse = 5*noeuds
394     
395      ! 1. mesh building
396      ! ----------------
397     
398      ! definition of specific information for a subdomain
399      !  narea           : subdomain number = processor number +1
400      !  ninterf         : neighbour subdomain number
401      !  nni             : interface point number
402      !  ndvois array    : neighbour subdomain list
403      !  maplistin array : node pointer at each interface
404      !  maplistin array : concatened list of interface nodes
405     
406      !  messag coding is necessary by interface type for avoid collision
407      !  if nperio == 1
408     
409      !  lint  array     : indoor interface list / type
410      !  lext  array     : outdoor interface list / type
411     
412      !  domain with jpniXjpnj subdomains
413     
414      CALL feti_inisub(nifmat,njfmat,nbondi,nbondj,nperio,   &
415          nbsw,nbnw,nbse,nbne,ninterf,ninterfc,nni,nnic)
416
417      CALL feti_creadr(malim,malimax,nim,3*ninterf ,mandvois ,'ndvois' )
418      CALL feti_creadr(malim,malimax,nim,3*ninterfc,mandvoisc,'ndvoisc')
419      CALL feti_creadr(malim,malimax,nim,ninterfc+1,maplistin,'plistin')
420      CALL feti_creadr(malim,malimax,nim,nnic      ,malistin ,'listin' )
421
422      ! pb with periodicity conditions : iiend, ijend
423
424      ijend = nlej
425      iiend = nlei
426      IF (jperio == 1) iiend = nlci - jpreci
427
428      CALL feti_subound(nifmat,njfmat,nldi,iiend,nldj,ijend,   &
429          narea,nbondi,nbondj,nperio,   &
430          ninterf,ninterfc,   &
431          nowe,noea,noso,nono,   &
432          nbsw,nbnw,nbse,nbne,   &
433          npsw,npnw,npse,npne,   &
434          mfet(mandvois),mfet(mandvoisc),   &
435          mfet(maplistin),nnic,mfet(malistin) )
436
437   END SUBROUTINE fetstr
438
439
440   SUBROUTINE fetmat
441      !!---------------------------------------------------------------------
442      !!                   ***  ROUTINE fetmat  ***
443      !!           
444      !! ** Purpose :   Construction of the matrix of the barotropic stream
445      !!      function system.
446      !!      Finite Elements Tearing & Interconnecting (FETI) approach
447      !!      Matrix treatment : Neumann condition, inverse computation
448      !!
449      !! ** Method :
450      !!
451      !! References :
452      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
453      !!      A domain decomposition solver to compute the barotropic
454      !!      component of an OGCM in the parallel processing field.
455      !!      Ocean Modelling, issue 105, december 94.
456      !!
457      !! History :
458      !!        !  98-02 (M. Guyon)  Original code
459      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
460      !!----------------------------------------------------------------------
461      !! * Modules used
462      USE lib_feti            ! feti librairy
463      !! * Local declarations
464      INTEGER ::   ji, jj, jk, jl
465      INTEGER ::   iimask(jpi,jpj)
466      INTEGER ::   iiend, ijend
467      REAL(wp) ::   zres, zres2, zdemi
468      !!---------------------------------------------------------------------
469
470      ! Matrix computation
471      ! ------------------
472
473      CALL feti_creadr(malxm,malxmax,nxm,nmorse,maan,'matrice a')
474
475      nnitot = nni
476
477      CALL mpp_sum( nnitot, 1, numit0ete )
478      CALL feti_creadr(malxm,malxmax,nxm,npe*npe,maae,'ae')
479
480      ! initialisation of the local barotropic matrix
481      ! local boundary conditions on the halo
482
483      CALL lbc_lnk( gcp(:,:,1), 'F', 1)
484      CALL lbc_lnk( gcp(:,:,2), 'F', 1) 
485      CALL lbc_lnk( gcp(:,:,3), 'F', 1) 
486      CALL lbc_lnk( gcp(:,:,4), 'F', 1) 
487      CALL lbc_lnk( gcdmat    , 'T', 1)
488
489      ! Neumann conditions
490      ! initialisation of the integer Neumann Mask
491
492      CALL feti_iclr(jpi*jpj,iimask)
493      DO jj = 1, jpj
494         DO ji = 1, jpi
495            iimask(ji,jj) = INT( gcdprc(ji,jj) )
496         END DO
497      END DO
498
499      ! regularization of the local matrix
500
501      DO jj = 1, jpj
502         DO ji = 1, jpi
503            gcdmat(ji,jj) = gcdmat(ji,jj) * gcdprc(ji,jj) + 1. - gcdprc(ji,jj)
504         END DO
505      END DO
506
507      DO jk = 1, 4
508         DO jj = 1, jpj
509            DO ji = 1, jpi
510               gcp(ji,jj,jk) = gcp(ji,jj,jk) * gcdprc(ji,jj)
511            END DO
512         END DO
513      END DO
514     
515      ! implementation of the west, east, north & south Neumann conditions
516
517      zdemi  = 0.5
518
519      ! pb with periodicity conditions : iiend, ijend
520
521      ijend = nlej
522      iiend = nlei
523      IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) iiend = nlci - jpreci
524
525      IF( nbondi == 2 .AND. (nperio /= 1 .OR. nperio /= 4 .OR. nperio == 6) ) THEN
526
527         ! with the periodicity : no east/west interface if nbondi = 2
528         ! and nperio != 1
529
530      ELSE 
531         ! west
532         IF( nbondi /= -1 ) THEN
533            DO jj = 1, jpj
534               IF( iimask(1,jj) /= 0 ) THEN
535                  gcp(1,jj,2) = 0.e0
536                  gcp(1,jj,1) = zdemi * gcp(1,jj,1)
537                  gcp(1,jj,4) = zdemi * gcp(1,jj,4)
538               ENDIF
539            END DO
540            DO jj = 1, jpj
541               IF( iimask(1,jj) /= 0 ) THEN
542                  gcdmat(1,jj) = - ( gcp(1,jj,1) + gcp(1,jj,2) + gcp(1,jj,3) + gcp(1,jj,4) )
543               ENDIF
544            END DO
545         ENDIF
546         ! east
547         IF( nbondi /= 1 ) THEN
548            DO jj = 1, jpj
549               IF( iimask(iiend,jj) /= 0 ) THEN
550                  gcp(iiend,jj,3) = 0.e0
551                  gcp(iiend,jj,1) = zdemi * gcp(iiend,jj,1)
552                  gcp(iiend,jj,4) = zdemi * gcp(iiend,jj,4)
553               ENDIF
554            END DO
555            DO jj = 1, jpj
556               IF( iimask(iiend,jj) /= 0 ) THEN
557                  gcdmat(iiend,jj) = - ( gcp(iiend,jj,1) + gcp(iiend,jj,2)   &
558                                       + gcp(iiend,jj,3) + gcp(iiend,jj,4) )
559               ENDIF
560            END DO
561         ENDIF
562      ENDIF
563
564      ! south
565      IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN
566         DO ji = 1, jpi
567            IF( iimask(ji,1) /= 0 ) THEN
568               gcp(ji,1,1) = 0.e0
569               gcp(ji,1,2) = zdemi * gcp(ji,1,2)
570               gcp(ji,1,3) = zdemi * gcp(ji,1,3)
571            ENDIF
572         END DO
573         DO ji = 1, jpi
574            IF( iimask(ji,1) /= 0 ) THEN
575               gcdmat(ji,1) = - ( gcp(ji,1,1) + gcp(ji,1,2) + gcp(ji,1,3) + gcp(ji,1,4) )
576            ENDIF
577         END DO
578      ENDIF
579     
580      ! north
581      IF( nbondj /= 1 .AND. nbondj /= 2 ) THEN
582         DO ji = 1, jpi
583            IF( iimask(ji,ijend) /= 0 ) THEN
584               gcp(ji,ijend,4) = 0.e0
585               gcp(ji,ijend,2) = zdemi * gcp(ji,ijend,2) 
586               gcp(ji,ijend,3) = zdemi * gcp(ji,ijend,3)
587            ENDIF
588         END DO
589         DO ji = 1, jpi
590            IF( iimask(ji,ijend) /= 0 ) THEN
591               gcdmat(ji,ijend) = - ( gcp(ji,ijend,1) + gcp(ji,ijend,2)   &
592                                    + gcp(ji,ijend,3) + gcp(ji,ijend,4) )
593            ENDIF
594         END DO
595      ENDIF
596
597      ! matrix terms are  saved in FETI solver arrays
598      CALL feti_vmov(noeuds,gcp(1,1,1),wfeti(maan))
599      CALL feti_vmov(noeuds,gcp(1,1,2),wfeti(maan+noeuds))
600      CALL feti_vmov(noeuds,gcdmat,wfeti(maan+2*noeuds))
601      CALL feti_vmov(noeuds,gcp(1,1,3),wfeti(maan+3*noeuds))
602      CALL feti_vmov(noeuds,gcp(1,1,4),wfeti(maan+4*noeuds))
603
604      ! construction of Dirichlet liberty degrees array
605      CALL feti_subdir(nifmat,njfmat,noeuds,ndir,iimask)
606      CALL feti_creadr(malim,malimax,nim,ndir,malisdir,'lisdir')
607      CALL feti_listdir(jpi,jpj,iimask,ndir,mfet(malisdir))
608
609      ! stop onto  matrix term for Dirichlet conditions
610      CALL feti_blomat(nifmat+1,njfmat+1,wfeti(maan),ndir,mfet(malisdir))
611
612      ! reservation of factorized diagonal blocs and temporary array for
613      ! factorization
614      npblo = (njfmat+1) * (nifmat+1) * (nifmat+1)
615      ndimax = nifmat+1
616
617      CALL feti_creadr(malxm,malxmax,nxm,npblo,mablo,'blo')
618      CALL feti_creadr(malxm,malxmax,nxm,noeuds,madia,'dia')
619      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mav,'v')
620      CALL feti_creadr(malxm,malxmax,nxm,ndimax*ndimax,mautil,'util')
621
622      ! stoping the rigid modes
623
624      ! the number of rigid modes =< Max [dim(Ker(Ep))]
625      !                                p=1,Np
626
627      CALL feti_creadr(malim,malimax,nim,ndkerep,malisblo,'lisblo')
628
629      ! Matrix factorization
630
631      CALL feti_front(noeuds,nifmat+1,njfmat+1,wfeti(maan),npblo,   &
632          wfeti(mablo),wfeti(madia),   &
633          wfeti(mautil),wfeti(mav),ndlblo,mfet(malisblo),ndkerep)
634      CALL feti_prext(noeuds,wfeti(madia))
635
636      ! virtual dealloc => we have to see for a light f90 version
637      ! the super structure is removed to clean the coarse grid
638      ! solver structure
639
640      malxm = madia
641      CALL feti_vclr(noeuds,wfeti(madia))
642      CALL feti_vclr(noeuds,wfeti(mav))
643      CALL feti_vclr(ndimax*ndimax,wfeti(mautil))
644
645      ! ndlblo is the dimension of the local nullspace .=<. the size of the
646      ! memory of the superstructure associated to the nullspace : ndkerep
647      ! ndkerep is introduced to avoid messages "out of bounds" when memory
648      ! is checked
649
650      ! copy matrix for Dirichlet condition
651
652      CALL feti_creadr(malxm,malxmax,nxm,noeuds,miax,'x')
653      CALL feti_creadr(malxm,malxmax,nxm,noeuds,may,'y')
654      CALL feti_creadr(malxm,malxmax,nxm,noeuds,maz,'z')
655
656      ! stoping the rigid modes
657
658      ! ndlblo is the dimension of the local nullspace .=<. the size of the
659      ! memory of the superstructure associated to the nullspace : ndkerep
660      ! ndkerep is introduced to avoid messages "out of bounds" when memory
661      ! is checked
662
663      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*noeuds,mansp,'nsp')
664      CALL feti_blomat1(nifmat+1,njfmat+1,wfeti(maan),ndlblo,   &
665          mfet(malisblo),wfeti(mansp))     
666
667      ! computation of operator kernel
668
669      CALL feti_nullsp(noeuds,nifmat+1,njfmat+1,npblo,wfeti(mablo),   &
670          wfeti(maan),ndlblo,mfet(malisblo),wfeti(mansp),   &
671          wfeti(maz))
672
673      ! test of the factorisation onto each sub domain
674
675      CALL feti_init(noeuds,wfeti(may))
676      CALL feti_blodir(noeuds,wfeti(may),ndir,mfet(malisdir))
677      CALL feti_blodir(noeuds,wfeti(may),ndlblo,mfet(malisblo))
678      CALL feti_vclr(noeuds,wfeti(miax))
679      CALL feti_resloc(noeuds,nifmat+1,njfmat+1,wfeti(maan),npblo,   &
680          wfeti(mablo),wfeti(may),wfeti(miax),wfeti(maz)) 
681      CALL feti_proax(noeuds,nifmat+1,njfmat+1,wfeti(maan),wfeti(miax),   &
682          wfeti(maz))
683      CALL feti_blodir(noeuds,wfeti(maz),ndlblo,mfet(malisblo))
684      CALL feti_vsub(noeuds,wfeti(may),wfeti(maz),wfeti(maz))
685
686      zres2 = 0.e0
687      DO jl = 1, noeuds
688         zres2 = zres2 + wfeti(may+jl-1) * wfeti(may+jl-1)
689      END DO
690      CALL mpp_sum(zres2,1,zres)
691
692      res2 = 0.e0
693      DO jl = 1, noeuds
694         res2 = res2 + wfeti(maz+jl-1) * wfeti(maz+jl-1)
695      END DO
696      res2 = res2 / zres2
697      CALL mpp_sum(res2,1,zres)
698
699      res2 = SQRT(res2)
700      IF(lwp) WRITE(numout,*) 'global residu : sqrt((Ax-b,Ax-b)/(b.b)) =', res2
701
702      IF( res2 > (eps/100.) ) THEN
703         IF(lwp) WRITE (numout,*) 'eps is :',eps
704         IF(lwp) WRITE (numout,*) 'factorized matrix precision :',res2
705         STOP
706      ENDIF
707
708   END SUBROUTINE fetmat
709
710
711   SUBROUTINE fetsch
712      !!---------------------------------------------------------------------
713      !!                  ***  ROUTINE fetsch  ***
714      !!       
715      !! ** Purpose :
716      !!     Construction of the matrix of the barotropic stream function
717      !!     system.
718      !!     Finite Elements Tearing & Interconnecting (FETI) approach
719      !!     Data framework for the Schur Dual solve
720      !!
721      !! ** Method :
722      !!
723      !! References :
724      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
725      !!      A domain decomposition solver to compute the barotropic
726      !!      component of an OGCM in the parallel processing field.
727      !!      Ocean Modelling, issue 105, december 94.
728      !!
729      !! History :
730      !!        !  98-02 (M. Guyon)  Original code
731      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
732      !!----------------------------------------------------------------------
733      !! * Modules used
734      USE lib_feti            ! feti librairy
735      !! * Local declarations
736      !!---------------------------------------------------------------------
737
738      !  computing weights for the conform construction
739
740      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mapoids ,'poids' )
741      CALL feti_creadr(malxm,malxmax,nxm,nnic  ,mabufin ,'bufin' )
742      CALL feti_creadr(malxm,malxmax,nxm,nnic  ,mabufout,'bufout')
743
744!!    CALL feti_poids(ninterfc,mfet(mandvoisc),mfet(maplistin),nnic,   &
745!!        mfet(malistin),narea,noeuds,wfeti(mapoids),wfeti(mabufin),   &
746!!        wfeti(mabufout) )
747      CALL feti_poids(ninterfc,                                nnic,   &
748          mfet(malistin),      noeuds,wfeti(mapoids)                )
749
750
751      ! Schur dual arrays
752 
753      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mabitw,'bitw') 
754      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mautilu,'utilu') 
755      CALL feti_creadr(malxm,malxmax,nxm,nni,malambda,'lambda') 
756      CALL feti_creadr(malxm,malxmax,nxm,nni,mag,'g') 
757      CALL feti_creadr(malxm,malxmax,nxm,nni,mapg,'pg') 
758      CALL feti_creadr(malxm,malxmax,nxm,nni,mamg,'mg') 
759      CALL feti_creadr(malxm,malxmax,nxm,nni,maw,'w') 
760      CALL feti_creadr(malxm,malxmax,nxm,nni,madw,'dw')
761
762      !  coarse grid solver dimension and arrays
763
764      nitmaxete = ndlblo
765      CALL  mpp_sum(nitmaxete,1,numit0ete)
766
767      nitmaxete = nitmaxete + 1
768      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maxnul,'xnul')
769      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maynul,'ynul')
770      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maeteg,'eteg')
771      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maeteag,'eteag')
772      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*nitmaxete,maeted,'eted')
773      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*nitmaxete,maetead,'etead')
774      CALL feti_creadr(malxm,malxmax,nxm,nitmaxete,maeteadd,'eteadd')
775      CALL feti_creadr(malxm,malxmax,nxm,nitmaxete,maetegamm,'etegamm')
776      CALL feti_creadr(malxm,malxmax,nxm,nni,maetew,'etew') 
777      CALL feti_creadr(malxm,malxmax,nxm,noeuds,maetev,'etev') 
778
779      ! construction of semi interface arrays
780
781      CALL feti_creadr(malim,malimax,nim,ninterf+1,maplistih,'plistih')
782!!    CALL feti_halfint(ninterf,mfet(mandvois),mfet(maplistin),nni,   &
783!!        mfet(maplistih),nnih,narea)
784      CALL feti_halfint(ninterf,mfet(mandvois),mfet(maplistin),       & 
785          mfet(maplistih),nnih      )
786
787      CALL feti_creadr(malxm,malxmax,nxm,nnih,magh,'gh') 
788
789      ! Schur Dual Method
790
791      nmaxd = nnitot / 2
792
793      !  computation of the remain array for descent directions
794
795      nmaxd = min0(nmaxd,(nxm-nitmaxete-malxm)/(2*nnih+3))
796      CALL mpp_min(nmaxd,1,numit0ete)
797
798      nitmax = nnitot/2
799      epsilo = eps
800      ntest = 0
801
802      ! Krylov space construction
803
804      CALL feti_creadr(malxm,malxmax,nxm,nnih*nmaxd,mawj,'wj') 
805      CALL feti_creadr(malxm,malxmax,nxm,nnih*nmaxd,madwj,'dwj') 
806      CALL feti_creadr(malxm,malxmax,nxm,nmaxd,madwwj,'dwwj') 
807      CALL feti_creadr(malxm,malxmax,nxm,nmaxd,magamm,'gamm') 
808      CALL feti_creadr(malxm,malxmax,nxm,max0(nmaxd,nitmaxete),mawork,'work')
809      mjj0 = 0 
810      numit0ete = 0 
811
812   END SUBROUTINE fetsch
813
814#else
815   SUBROUTINE fetstr                 ! Empty routine
816   END SUBROUTINE fetstr
817   SUBROUTINE fetmat                 ! Empty routine
818   END SUBROUTINE fetmat
819   SUBROUTINE fetsch                 ! Empty routine
820   END SUBROUTINE fetsch
821#endif
822
823   !!======================================================================
824END MODULE solmat
Note: See TracBrowser for help on using the repository browser.