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

Last change on this file since 359 was 359, checked in by opalod, 18 years ago

nemo_v1_update_033 : RB + CT : Add new surface pressure gradient algorithms
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.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
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_flt ) THEN   !bug missing lk_dynspg_flt_atsk
202#if ! defined key_dynspg_flt
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 ) .OR. ( nsolv == 4 ) )  gccd(:,:) = sor * gcp(:,:,2)
235
236         IF( nsolv == 4 ) THEN
237            CALL lbc_lnk_e( gcp   (:,:,1), c_solver_pt, 1. )   ! lateral boundary conditions
238            CALL lbc_lnk_e( gcp   (:,:,2), c_solver_pt, 1. )   ! lateral boundary conditions
239            CALL lbc_lnk_e( gcp   (:,:,3), c_solver_pt, 1. )   ! lateral boundary conditions
240            CALL lbc_lnk_e( gcp   (:,:,4), c_solver_pt, 1. )   ! lateral boundary conditions
241            CALL lbc_lnk_e( gcdprc(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions
242            CALL lbc_lnk_e( gcdmat(:,:)  , c_solver_pt, 1. )   ! lateral boundary conditions         
243            IF( npolj /= 0 ) CALL sol_exd( gcp , c_solver_pt ) ! switch northernelements
244         END IF
245
246      ELSE
247         
248         ! FETI method
249         ! if feti solver : gcdprc is a mask for the non-overlapping
250         !   data structuring
251         
252         DO jj = 1, jpj
253            DO ji = 1, jpi
254               IF( bmask(ji,jj) /= 0.e0 ) THEN
255                  gcdprc(ji,jj) = 1.e0
256               ELSE
257                  gcdprc(ji,jj) = 0.e0
258               ENDIF
259            END DO
260         END DO
261         
262         ! so "common" line & "common" column have to be !=0 except on global
263         !   domain boundaries
264         ! pbs with nbondi if nperio != 2 ?
265         !   ii = nldi-1
266         ! pb with nldi value if jperio==1 : nbondi modifyed at the end
267         !   of inimpp.F => pb
268         ! pb with periodicity conditions : iiend, ijend
269         
270         ijend = nlej
271         iiend = nlei
272         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 )   iiend = nlci - jpreci
273         ii = jpreci
274         
275         ! case number 1
276         
277         IF( nbondi /= -1 .AND. nbondi /= 2 ) THEN
278            DO jj = 1, ijend
279               IF( fmask(ii,jj,1) == 1. ) gcdprc(ii,jj) = 1.
280            END DO
281         ENDIF
282         
283         ! case number 2
284         
285         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
286            DO jj = 1, ijend
287               IF( fmask(ii,jj,1) == 1. ) gcdprc(ii,jj) = 1.
288            END DO
289         ENDIF
290         
291         !      ij = nldj-1
292         ! pb with nldi value if jperio==1 : nbondi modifyed at the end
293         !   of inimpp.F => pb, here homogeneisation...
294         
295         ij = jprecj
296         IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN
297            DO ji = 1, iiend
298               IF( fmask(ji,ij,1) == 1. ) gcdprc(ji,ij) = 1.
299            END DO
300         ENDIF
301      ENDIF
302     
303     
304      ! 4. Initialization the arrays used in pcg
305      ! ----------------------------------------
306      gcx  (:,:) = 0.e0
307      gcxb (:,:) = 0.e0
308      gcb  (:,:) = 0.e0
309      gcr  (:,:) = 0.e0
310      gcdes(:,:) = 0.e0
311      gccd (:,:) = 0.e0
312     
313      ! FETI method
314      IF( nsolv == 3 ) THEN
315         CALL fetmat       ! Matrix treatment : Neumann condition, inverse computation
316         CALL fetsch       ! data framework for the Schur Dual solver
317      ENDIF
318     
319   END SUBROUTINE sol_mat
320
321
322   SUBROUTINE sol_exd( pt3d, cd_type )
323      !!----------------------------------------------------------------------
324      !!                  ***  routine sol_exd  ***
325      !!                 
326      !! ** Purpose :   Reorder gcb coefficient on the extra outer  halo
327      !!                at north fold in case of T or F pivot
328      !!
329      !! ** Method  :   Perform a circular permutation of the coefficients on
330      !!                the total area strictly above the pivot point,
331      !!                and on the semi-row of the pivot point   
332      !!               
333      !! History :
334      !!   9.0  !  05-09  (R. Benshila)  original routine
335      !!----------------------------------------------------------------------
336      !! * Arguments
337      CHARACTER(len=1) , INTENT( in ) ::   &
338         cd_type       ! define the nature of pt2d array grid-points
339         !             !  = T , U , V , F , W
340         !             !  = S : T-point, north fold treatment
341         !             !  = G : F-point, north fold treatment
342         !             !  = I : sea-ice velocity at F-point with index shift
343      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4), INTENT( inout ) ::   &
344         pt3d          ! 2D array on which the boundary condition is applied
345
346      !! * Local variables
347      INTEGER  ::   ji, jk      ! dummy loop indices
348      INTEGER  ::   iloc                ! temporary integers
349      REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj,4) ::   &
350         ztab          ! 2D array on which the boundary condition is applied
351      !!----------------------------------------------------------------------
352
353      ztab = pt3d
354
355      ! north fold treatment
356      ! -----------------------
357 
358      SELECT CASE ( npolj )
359         
360         CASE ( 3 , 4 )   !  T pivot
361         iloc = jpiglo/2 +1 
362           
363            SELECT CASE ( cd_type )
364 
365            CASE ( 'T', 'S', 'U', 'W' )
366               DO jk =1, 4
367                  DO ji = 1-jpr2di, nlci+jpr2di
368                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
369                  ENDDO
370               ENDDO
371
372              DO jk =1, 4
373                  DO ji = nlci+jpr2di, 1-jpr2di,  -1
374                     IF( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) &
375                       & .OR. ( mi0(iloc) == jpi+1 ) ) EXIT
376                     pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
377                  ENDDO
378               ENDDO
379
380            CASE ( 'F' ,'G' , 'I', 'V' )
381               DO jk =1, 4
382                  DO ji = 1-jpr2di, nlci+jpr2di
383                     pt3d(ji,nlcj-1:nlcj+jpr2dj,jk) = ztab(ji,nlcj-1:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
384                  ENDDO
385               ENDDO
386
387            END SELECT   ! cd_type
388 
389         CASE ( 5 , 6 )                 ! F pivot
390          iloc=jpiglo/2
391
392            SELECT CASE (cd_type )
393
394            CASE ( 'T'  ,'S', 'U', 'W')
395               DO jk =1, 4
396                  DO ji = 1-jpr2di, nlci+jpr2di
397                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
398                  ENDDO
399               ENDDO
400
401            CASE ( 'F' ,'G' , 'I', 'V' )
402               DO jk =1, 4
403                  DO ji = 1-jpr2di, nlci+jpr2di
404                     pt3d(ji,nlcj:nlcj+jpr2dj,jk) = ztab(ji,nlcj:nlcj+jpr2dj,jk+3-2*MOD(jk+3,4))           
405                  ENDDO
406               ENDDO
407               DO jk =1, 4
408                  DO ji = nlci+jpr2di, 1-jpr2di,  -1
409                    IF ( ( ji .LT. mi0(iloc) .AND. mi0(iloc) /= 1 ) &
410                       & .OR. ( mi0(iloc) == jpi+1 ) ) EXIT
411                    pt3d(ji,nlcj-1,jk) = ztab(ji,nlcj-1,jk+3-2*MOD(jk+3,4))
412                  ENDDO
413               ENDDO
414
415            END SELECT   ! cd_type
416
417         END SELECT   ! npolj
418 
419   END SUBROUTINE sol_exd
420
421#if defined key_feti
422
423   SUBROUTINE fetstr
424      !!---------------------------------------------------------------------
425      !!                  ***  ROUTINE fetstr  ***
426      !!               
427      !! ** Purpose :   Construction of the matrix of the barotropic stream
428      !!      function system.
429      !!      Finite Elements Tearing & Interconnecting (FETI) approach
430      !!      Memory allocation and interface definition for the solver
431      !!
432      !! ** Method :
433      !!
434      !! References :
435      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
436      !!      A domain decomposition solver to compute the barotropic
437      !!      component of an OGCM in the parallel processing field.
438      !!      Ocean Modelling, issue 105, december 94.
439      !!
440      !! History :
441      !!        !  98-02 (M. Guyon)  Original code
442      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
443      !!----------------------------------------------------------------------
444      !! * Modules used
445      USE lib_feti            ! feti librairy
446      !! * Local declarations
447      INTEGER ::   iiend, ijend, iperio   ! temporary integers
448      !!---------------------------------------------------------------------
449     
450     
451      ! Preconditioning technics of the Dual Schur Operator
452      ! <= definition of the Coarse Grid solver
453      ! <= dimension of the nullspace of the local operators
454      ! <= Neumann boundaries conditions
455     
456      ! 0. Initializations
457      ! ------------------
458     
459      ndkerep = 1
460
461      ! initialization of the superstructures management
462
463      malxm = 1
464      malim = 1
465
466      ! memory space for the pcpg associated with the FETI dual formulation
467      ! ndkerep is associated to the list of rigid modes,
468      ! ndkerep == 1 because the Dual Operator
469      ! is a first order operator due to SPG elliptic Operator is a
470      ! second order operator
471
472      nim = 50
473      nim = nim + ndkerep
474      nim = nim + 2*jpi + 2*jpj
475      nim = nim + jpi*jpj
476
477      nxm = 33
478      nxm = nxm + 4*jpnij
479      nxm = nxm + 19*(jpi+jpj)
480      nxm = nxm + 13*jpi*jpj
481      nxm = nxm + jpi*jpi*jpj
482
483      ! krylov space memory
484
485      iperio = 0
486      IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6) iperio = 1
487      nxm = nxm + 3*(jpnij-jpni)*jpi
488      nxm = nxm + 3*(jpnij-jpnj+iperio)*jpj
489      nxm = nxm + 2*(jpi+jpj)*(jpnij-jpni)*jpi
490      nxm = nxm + 2*(jpi+jpj)*(jpnij-jpnj+iperio)*jpj
491
492      ! Resolution with the Schur dual Method ( frontal and local solver by
493      ! blocks
494      ! Case with a local symetrical matrix
495      ! The local matrix is stored in a multi-column form
496      ! The total number of nodes for this subdomain is named "noeuds"
497
498      noeuds = jpi*jpj
499      nifmat = jpi-1
500      njfmat = jpj-1
501      nelem = nifmat*njfmat
502      npe = 4
503      nmorse = 5*noeuds
504     
505      ! 1. mesh building
506      ! ----------------
507     
508      ! definition of specific information for a subdomain
509      !  narea           : subdomain number = processor number +1
510      !  ninterf         : neighbour subdomain number
511      !  nni             : interface point number
512      !  ndvois array    : neighbour subdomain list
513      !  maplistin array : node pointer at each interface
514      !  maplistin array : concatened list of interface nodes
515     
516      !  messag coding is necessary by interface type for avoid collision
517      !  if nperio == 1
518     
519      !  lint  array     : indoor interface list / type
520      !  lext  array     : outdoor interface list / type
521     
522      !  domain with jpniXjpnj subdomains
523     
524      CALL feti_inisub(nifmat,njfmat,nbondi,nbondj,nperio,   &
525          nbsw,nbnw,nbse,nbne,ninterf,ninterfc,nni,nnic)
526
527      CALL feti_creadr(malim,malimax,nim,3*ninterf ,mandvois ,'ndvois' )
528      CALL feti_creadr(malim,malimax,nim,3*ninterfc,mandvoisc,'ndvoisc')
529      CALL feti_creadr(malim,malimax,nim,ninterfc+1,maplistin,'plistin')
530      CALL feti_creadr(malim,malimax,nim,nnic      ,malistin ,'listin' )
531
532      ! pb with periodicity conditions : iiend, ijend
533
534      ijend = nlej
535      iiend = nlei
536      IF (jperio == 1) iiend = nlci - jpreci
537
538      CALL feti_subound(nifmat,njfmat,nldi,iiend,nldj,ijend,   &
539          narea,nbondi,nbondj,nperio,   &
540          ninterf,ninterfc,   &
541          nowe,noea,noso,nono,   &
542          nbsw,nbnw,nbse,nbne,   &
543          npsw,npnw,npse,npne,   &
544          mfet(mandvois),mfet(mandvoisc),   &
545          mfet(maplistin),nnic,mfet(malistin) )
546
547   END SUBROUTINE fetstr
548
549
550   SUBROUTINE fetmat
551      !!---------------------------------------------------------------------
552      !!                   ***  ROUTINE fetmat  ***
553      !!           
554      !! ** Purpose :   Construction of the matrix of the barotropic stream
555      !!      function system.
556      !!      Finite Elements Tearing & Interconnecting (FETI) approach
557      !!      Matrix treatment : Neumann condition, inverse computation
558      !!
559      !! ** Method :
560      !!
561      !! References :
562      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
563      !!      A domain decomposition solver to compute the barotropic
564      !!      component of an OGCM in the parallel processing field.
565      !!      Ocean Modelling, issue 105, december 94.
566      !!
567      !! History :
568      !!        !  98-02 (M. Guyon)  Original code
569      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
570      !!----------------------------------------------------------------------
571      !! * Modules used
572      USE lib_feti            ! feti librairy
573      !! * Local declarations
574      INTEGER ::   ji, jj, jk, jl
575      INTEGER ::   iimask(jpi,jpj)
576      INTEGER ::   iiend, ijend
577      REAL(wp) ::   zres, zres2, zdemi
578      !!---------------------------------------------------------------------
579
580      ! Matrix computation
581      ! ------------------
582
583      CALL feti_creadr(malxm,malxmax,nxm,nmorse,maan,'matrice a')
584
585      nnitot = nni
586
587      CALL mpp_sum( nnitot, 1, numit0ete )
588      CALL feti_creadr(malxm,malxmax,nxm,npe*npe,maae,'ae')
589
590      ! initialisation of the local barotropic matrix
591      ! local boundary conditions on the halo
592
593      CALL lbc_lnk( gcp(:,:,1), 'F', 1)
594      CALL lbc_lnk( gcp(:,:,2), 'F', 1) 
595      CALL lbc_lnk( gcp(:,:,3), 'F', 1) 
596      CALL lbc_lnk( gcp(:,:,4), 'F', 1) 
597      CALL lbc_lnk( gcdmat    , 'T', 1)
598
599      ! Neumann conditions
600      ! initialisation of the integer Neumann Mask
601
602      CALL feti_iclr(jpi*jpj,iimask)
603      DO jj = 1, jpj
604         DO ji = 1, jpi
605            iimask(ji,jj) = INT( gcdprc(ji,jj) )
606         END DO
607      END DO
608
609      ! regularization of the local matrix
610
611      DO jj = 1, jpj
612         DO ji = 1, jpi
613            gcdmat(ji,jj) = gcdmat(ji,jj) * gcdprc(ji,jj) + 1. - gcdprc(ji,jj)
614         END DO
615      END DO
616
617      DO jk = 1, 4
618         DO jj = 1, jpj
619            DO ji = 1, jpi
620               gcp(ji,jj,jk) = gcp(ji,jj,jk) * gcdprc(ji,jj)
621            END DO
622         END DO
623      END DO
624     
625      ! implementation of the west, east, north & south Neumann conditions
626
627      zdemi  = 0.5
628
629      ! pb with periodicity conditions : iiend, ijend
630
631      ijend = nlej
632      iiend = nlei
633      IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) iiend = nlci - jpreci
634
635      IF( nbondi == 2 .AND. (nperio /= 1 .OR. nperio /= 4 .OR. nperio == 6) ) THEN
636
637         ! with the periodicity : no east/west interface if nbondi = 2
638         ! and nperio != 1
639
640      ELSE 
641         ! west
642         IF( nbondi /= -1 ) THEN
643            DO jj = 1, jpj
644               IF( iimask(1,jj) /= 0 ) THEN
645                  gcp(1,jj,2) = 0.e0
646                  gcp(1,jj,1) = zdemi * gcp(1,jj,1)
647                  gcp(1,jj,4) = zdemi * gcp(1,jj,4)
648               ENDIF
649            END DO
650            DO jj = 1, jpj
651               IF( iimask(1,jj) /= 0 ) THEN
652                  gcdmat(1,jj) = - ( gcp(1,jj,1) + gcp(1,jj,2) + gcp(1,jj,3) + gcp(1,jj,4) )
653               ENDIF
654            END DO
655         ENDIF
656         ! east
657         IF( nbondi /= 1 ) THEN
658            DO jj = 1, jpj
659               IF( iimask(iiend,jj) /= 0 ) THEN
660                  gcp(iiend,jj,3) = 0.e0
661                  gcp(iiend,jj,1) = zdemi * gcp(iiend,jj,1)
662                  gcp(iiend,jj,4) = zdemi * gcp(iiend,jj,4)
663               ENDIF
664            END DO
665            DO jj = 1, jpj
666               IF( iimask(iiend,jj) /= 0 ) THEN
667                  gcdmat(iiend,jj) = - ( gcp(iiend,jj,1) + gcp(iiend,jj,2)   &
668                                       + gcp(iiend,jj,3) + gcp(iiend,jj,4) )
669               ENDIF
670            END DO
671         ENDIF
672      ENDIF
673
674      ! south
675      IF( nbondj /= -1 .AND. nbondj /= 2 ) THEN
676         DO ji = 1, jpi
677            IF( iimask(ji,1) /= 0 ) THEN
678               gcp(ji,1,1) = 0.e0
679               gcp(ji,1,2) = zdemi * gcp(ji,1,2)
680               gcp(ji,1,3) = zdemi * gcp(ji,1,3)
681            ENDIF
682         END DO
683         DO ji = 1, jpi
684            IF( iimask(ji,1) /= 0 ) THEN
685               gcdmat(ji,1) = - ( gcp(ji,1,1) + gcp(ji,1,2) + gcp(ji,1,3) + gcp(ji,1,4) )
686            ENDIF
687         END DO
688      ENDIF
689     
690      ! north
691      IF( nbondj /= 1 .AND. nbondj /= 2 ) THEN
692         DO ji = 1, jpi
693            IF( iimask(ji,ijend) /= 0 ) THEN
694               gcp(ji,ijend,4) = 0.e0
695               gcp(ji,ijend,2) = zdemi * gcp(ji,ijend,2) 
696               gcp(ji,ijend,3) = zdemi * gcp(ji,ijend,3)
697            ENDIF
698         END DO
699         DO ji = 1, jpi
700            IF( iimask(ji,ijend) /= 0 ) THEN
701               gcdmat(ji,ijend) = - ( gcp(ji,ijend,1) + gcp(ji,ijend,2)   &
702                                    + gcp(ji,ijend,3) + gcp(ji,ijend,4) )
703            ENDIF
704         END DO
705      ENDIF
706
707      ! matrix terms are  saved in FETI solver arrays
708      CALL feti_vmov(noeuds,gcp(1,1,1),wfeti(maan))
709      CALL feti_vmov(noeuds,gcp(1,1,2),wfeti(maan+noeuds))
710      CALL feti_vmov(noeuds,gcdmat,wfeti(maan+2*noeuds))
711      CALL feti_vmov(noeuds,gcp(1,1,3),wfeti(maan+3*noeuds))
712      CALL feti_vmov(noeuds,gcp(1,1,4),wfeti(maan+4*noeuds))
713
714      ! construction of Dirichlet liberty degrees array
715      CALL feti_subdir(nifmat,njfmat,noeuds,ndir,iimask)
716      CALL feti_creadr(malim,malimax,nim,ndir,malisdir,'lisdir')
717      CALL feti_listdir(jpi,jpj,iimask,ndir,mfet(malisdir))
718
719      ! stop onto  matrix term for Dirichlet conditions
720      CALL feti_blomat(nifmat+1,njfmat+1,wfeti(maan),ndir,mfet(malisdir))
721
722      ! reservation of factorized diagonal blocs and temporary array for
723      ! factorization
724      npblo = (njfmat+1) * (nifmat+1) * (nifmat+1)
725      ndimax = nifmat+1
726
727      CALL feti_creadr(malxm,malxmax,nxm,npblo,mablo,'blo')
728      CALL feti_creadr(malxm,malxmax,nxm,noeuds,madia,'dia')
729      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mav,'v')
730      CALL feti_creadr(malxm,malxmax,nxm,ndimax*ndimax,mautil,'util')
731
732      ! stoping the rigid modes
733
734      ! the number of rigid modes =< Max [dim(Ker(Ep))]
735      !                                p=1,Np
736
737      CALL feti_creadr(malim,malimax,nim,ndkerep,malisblo,'lisblo')
738
739      ! Matrix factorization
740
741      CALL feti_front(noeuds,nifmat+1,njfmat+1,wfeti(maan),npblo,   &
742          wfeti(mablo),wfeti(madia),   &
743          wfeti(mautil),wfeti(mav),ndlblo,mfet(malisblo),ndkerep)
744      CALL feti_prext(noeuds,wfeti(madia))
745
746      ! virtual dealloc => we have to see for a light f90 version
747      ! the super structure is removed to clean the coarse grid
748      ! solver structure
749
750      malxm = madia
751      CALL feti_vclr(noeuds,wfeti(madia))
752      CALL feti_vclr(noeuds,wfeti(mav))
753      CALL feti_vclr(ndimax*ndimax,wfeti(mautil))
754
755      ! ndlblo is the dimension of the local nullspace .=<. the size of the
756      ! memory of the superstructure associated to the nullspace : ndkerep
757      ! ndkerep is introduced to avoid messages "out of bounds" when memory
758      ! is checked
759
760      ! copy matrix for Dirichlet condition
761
762      CALL feti_creadr(malxm,malxmax,nxm,noeuds,miax,'x')
763      CALL feti_creadr(malxm,malxmax,nxm,noeuds,may,'y')
764      CALL feti_creadr(malxm,malxmax,nxm,noeuds,maz,'z')
765
766      ! stoping the rigid modes
767
768      ! ndlblo is the dimension of the local nullspace .=<. the size of the
769      ! memory of the superstructure associated to the nullspace : ndkerep
770      ! ndkerep is introduced to avoid messages "out of bounds" when memory
771      ! is checked
772
773      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*noeuds,mansp,'nsp')
774      CALL feti_blomat1(nifmat+1,njfmat+1,wfeti(maan),ndlblo,   &
775          mfet(malisblo),wfeti(mansp))     
776
777      ! computation of operator kernel
778
779      CALL feti_nullsp(noeuds,nifmat+1,njfmat+1,npblo,wfeti(mablo),   &
780          wfeti(maan),ndlblo,mfet(malisblo),wfeti(mansp),   &
781          wfeti(maz))
782
783      ! test of the factorisation onto each sub domain
784
785      CALL feti_init(noeuds,wfeti(may))
786      CALL feti_blodir(noeuds,wfeti(may),ndir,mfet(malisdir))
787      CALL feti_blodir(noeuds,wfeti(may),ndlblo,mfet(malisblo))
788      CALL feti_vclr(noeuds,wfeti(miax))
789      CALL feti_resloc(noeuds,nifmat+1,njfmat+1,wfeti(maan),npblo,   &
790          wfeti(mablo),wfeti(may),wfeti(miax),wfeti(maz)) 
791      CALL feti_proax(noeuds,nifmat+1,njfmat+1,wfeti(maan),wfeti(miax),   &
792          wfeti(maz))
793      CALL feti_blodir(noeuds,wfeti(maz),ndlblo,mfet(malisblo))
794      CALL feti_vsub(noeuds,wfeti(may),wfeti(maz),wfeti(maz))
795
796      zres2 = 0.e0
797      DO jl = 1, noeuds
798         zres2 = zres2 + wfeti(may+jl-1) * wfeti(may+jl-1)
799      END DO
800      CALL mpp_sum(zres2,1,zres)
801
802      res2 = 0.e0
803      DO jl = 1, noeuds
804         res2 = res2 + wfeti(maz+jl-1) * wfeti(maz+jl-1)
805      END DO
806      res2 = res2 / zres2
807      CALL mpp_sum(res2,1,zres)
808
809      res2 = SQRT(res2)
810      IF(lwp) WRITE(numout,*) 'global residu : sqrt((Ax-b,Ax-b)/(b.b)) =', res2
811
812      IF( res2 > (eps/100.) ) THEN
813         IF(lwp) WRITE (numout,*) 'eps is :',eps
814         IF(lwp) WRITE (numout,*) 'factorized matrix precision :',res2
815         STOP
816      ENDIF
817
818   END SUBROUTINE fetmat
819
820
821   SUBROUTINE fetsch
822      !!---------------------------------------------------------------------
823      !!                  ***  ROUTINE fetsch  ***
824      !!       
825      !! ** Purpose :
826      !!     Construction of the matrix of the barotropic stream function
827      !!     system.
828      !!     Finite Elements Tearing & Interconnecting (FETI) approach
829      !!     Data framework for the Schur Dual solve
830      !!
831      !! ** Method :
832      !!
833      !! References :
834      !!      Guyon, M, Roux, F-X, Chartier, M and Fraunie, P, 1994 :
835      !!      A domain decomposition solver to compute the barotropic
836      !!      component of an OGCM in the parallel processing field.
837      !!      Ocean Modelling, issue 105, december 94.
838      !!
839      !! History :
840      !!        !  98-02 (M. Guyon)  Original code
841      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
842      !!----------------------------------------------------------------------
843      !! * Modules used
844      USE lib_feti            ! feti librairy
845      !! * Local declarations
846      !!---------------------------------------------------------------------
847
848      !  computing weights for the conform construction
849
850      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mapoids ,'poids' )
851      CALL feti_creadr(malxm,malxmax,nxm,nnic  ,mabufin ,'bufin' )
852      CALL feti_creadr(malxm,malxmax,nxm,nnic  ,mabufout,'bufout')
853
854!!    CALL feti_poids(ninterfc,mfet(mandvoisc),mfet(maplistin),nnic,   &
855!!        mfet(malistin),narea,noeuds,wfeti(mapoids),wfeti(mabufin),   &
856!!        wfeti(mabufout) )
857      CALL feti_poids(ninterfc,                                nnic,   &
858          mfet(malistin),      noeuds,wfeti(mapoids)                )
859
860
861      ! Schur dual arrays
862 
863      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mabitw,'bitw') 
864      CALL feti_creadr(malxm,malxmax,nxm,noeuds,mautilu,'utilu') 
865      CALL feti_creadr(malxm,malxmax,nxm,nni,malambda,'lambda') 
866      CALL feti_creadr(malxm,malxmax,nxm,nni,mag,'g') 
867      CALL feti_creadr(malxm,malxmax,nxm,nni,mapg,'pg') 
868      CALL feti_creadr(malxm,malxmax,nxm,nni,mamg,'mg') 
869      CALL feti_creadr(malxm,malxmax,nxm,nni,maw,'w') 
870      CALL feti_creadr(malxm,malxmax,nxm,nni,madw,'dw')
871
872      !  coarse grid solver dimension and arrays
873
874      nitmaxete = ndlblo
875      CALL  mpp_sum(nitmaxete,1,numit0ete)
876
877      nitmaxete = nitmaxete + 1
878      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maxnul,'xnul')
879      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maynul,'ynul')
880      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maeteg,'eteg')
881      CALL feti_creadr(malxm,malxmax,nxm,ndkerep,maeteag,'eteag')
882      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*nitmaxete,maeted,'eted')
883      CALL feti_creadr(malxm,malxmax,nxm,ndkerep*nitmaxete,maetead,'etead')
884      CALL feti_creadr(malxm,malxmax,nxm,nitmaxete,maeteadd,'eteadd')
885      CALL feti_creadr(malxm,malxmax,nxm,nitmaxete,maetegamm,'etegamm')
886      CALL feti_creadr(malxm,malxmax,nxm,nni,maetew,'etew') 
887      CALL feti_creadr(malxm,malxmax,nxm,noeuds,maetev,'etev') 
888
889      ! construction of semi interface arrays
890
891      CALL feti_creadr(malim,malimax,nim,ninterf+1,maplistih,'plistih')
892!!    CALL feti_halfint(ninterf,mfet(mandvois),mfet(maplistin),nni,   &
893!!        mfet(maplistih),nnih,narea)
894      CALL feti_halfint(ninterf,mfet(mandvois),mfet(maplistin),       & 
895          mfet(maplistih),nnih      )
896
897      CALL feti_creadr(malxm,malxmax,nxm,nnih,magh,'gh') 
898
899      ! Schur Dual Method
900
901      nmaxd = nnitot / 2
902
903      !  computation of the remain array for descent directions
904
905      nmaxd = min0(nmaxd,(nxm-nitmaxete-malxm)/(2*nnih+3))
906      CALL mpp_min(nmaxd,1,numit0ete)
907
908      nitmax = nnitot/2
909      epsilo = eps
910      ntest = 0
911
912      ! Krylov space construction
913
914      CALL feti_creadr(malxm,malxmax,nxm,nnih*nmaxd,mawj,'wj') 
915      CALL feti_creadr(malxm,malxmax,nxm,nnih*nmaxd,madwj,'dwj') 
916      CALL feti_creadr(malxm,malxmax,nxm,nmaxd,madwwj,'dwwj') 
917      CALL feti_creadr(malxm,malxmax,nxm,nmaxd,magamm,'gamm') 
918      CALL feti_creadr(malxm,malxmax,nxm,max0(nmaxd,nitmaxete),mawork,'work')
919      mjj0 = 0 
920      numit0ete = 0 
921
922   END SUBROUTINE fetsch
923
924#else
925   SUBROUTINE fetstr                 ! Empty routine
926   END SUBROUTINE fetstr
927   SUBROUTINE fetmat                 ! Empty routine
928   END SUBROUTINE fetmat
929   SUBROUTINE fetsch                 ! Empty routine
930   END SUBROUTINE fetsch
931#endif
932
933   !!======================================================================
934END MODULE solmat
Note: See TracBrowser for help on using the repository browser.