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

source: tags/nemo_v1_07/NEMO/OPA_SRC/SOL/solmat.F90 @ 7829

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

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