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

Last change on this file since 111 was 111, checked in by opalod, 20 years ago

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