source: trunk/NEMO/OPA_SRC/SOL/solisl.F90 @ 746

Last change on this file since 746 was 746, checked in by smasson, 14 years ago

implement ldstop in iom_varid, see ticket:21

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 44.7 KB
Line 
1MODULE solisl
2   !!==============================================================================
3   !!                       ***  MODULE  solisl  ***
4   !! Ocean island : specific treatment of island in rigid-lid case
5   !!==============================================================================
6#if defined key_islands
7   !!----------------------------------------------------------------------
8   !!   'key_islands' :                           islands in rigid-lid case
9   !!----------------------------------------------------------------------
10   !!   isl_dom     : locate islands in the domain
11   !!   isl_pri     : control print of island gridpoint position
12   !!   isl_pth     : Compute coeff. associated with path round each island
13   !!   isl_mat     : Compute the matrix associated with islands
14   !!   isl_bsf     : Compute the barotropic streamfunction of each island
15   !!   isl_dyn_spg : Update the barotropic streamfunction trend with the
16   !!                 Island contribution (call by dyn_spg)
17   !!   isl_stp_ctl : print island information (call by stp_ctl)
18   !!----------------------------------------------------------------------
19   !! * Modules used
20   USE oce             ! ocean dynamics and tracers variables
21   USE dom_oce         ! ocean space and time domain variables
22   USE in_out_manager  ! I/O manager
23   USE sol_oce         ! ocean solver
24   USE obc_oce         ! ocean open boundary condition
25   USE lib_mpp         ! distributed memory computing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Routine accessibility
31   PUBLIC isl_dom        ! routine called by solver_init
32   PUBLIC isl_mat        ! routine called by solver_init
33   PUBLIC isl_bsf        ! routine called by solver_init
34   PUBLIC isl_dyn_spg    ! routine called by dyn_spg
35   PUBLIC isl_stp_ctl    ! routine called by stp_ctl
36
37   !! * Shared module variables
38   LOGICAL, PUBLIC, PARAMETER ::   lk_isl = .TRUE.    !: 'key_islands' flag
39
40   !! * module variable
41   INTEGER :: numisl           ! logical unit for island file only used
42   !                           ! here during the initialization phase
43   INTEGER ::   &
44      nimlu, njmlu, nkmlu,  &  ! i-j-k-dimensions read
45      nlmlu, nmmlu, nnmlu      ! read islands number
46   INTEGER, DIMENSION(jpnisl,0:4,jpisl) ::   &
47      miisl, mjisl             ! position of island grid-points
48   INTEGER, DIMENSION(0:4,jpisl) ::   &
49      mnisl                    ! number of grid-points for each island
50 
51   REAL(wp) ::    &
52      replu                    ! read absolute precision
53   REAL(wp), DIMENSION(jpisl,jpisl) ::   &
54      aisl, aislm1             ! island matrix and its inverse
55   REAL(wp), DIMENSION(jpi,jpj,jpisl) ::   &
56      bsfisl                   ! barotropic streamfunction of island
57   REAL(wp), DIMENSION(jpi,jpj,2,jpisl) ::   &
58      acisl1, acisl2           ! coef. to compute circulations round islands
59   REAL(wp), DIMENSION(jpisl) ::   &
60      bisl,                 &  ! second member of island linear system
61      visl                     ! trend of island stream function
62   !!----------------------------------------------------------------------
63   !!   OPA 9.0 , LOCEAN-IPSL (2005)
64   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/SOL/solisl.F90,v 1.8 2007/06/05 10:40:06 opalod Exp $
65   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
66   !!----------------------------------------------------------------------
67
68CONTAINS
69
70   SUBROUTINE isl_dom
71      !!----------------------------------------------------------------------
72      !!                  ***  ROUTINE isl_dom  ***
73      !!                   
74      !! ** Purpose :   Locate island grid-points from mbathy array
75      !!
76      !! ** Method  :   The coordinates of ocean grid-points round an island
77      !!      are found from mbathy array read or computed in dommba.
78      !!      we first compute zwb, an ocean/land mask defined as follows:
79      !!              zwb(i,j)  =  0. over the main land and the ocean
80      !!                        = -n  over the nth island
81      !!      With the proper boundary conditions (defined by nperio)
82      !!        Note: IF i,j are the coordinates of an ocean grid-point west
83      !!      or east of a island, the corresponding coordinates miisl(ip,4,n)
84      !!      or mjisl(ip,2,n) are those of the western or southern side of
85      !!      the island (i.e. i+1 ou j+1, respectively)
86      !!
87      !! ** Action  :   compute mnisl, miisl, mjisl arrays defined as:
88      !!      mnisl(i,n) : nb of grid-points along (i=0), north (i=1), north
89      !!      (i=1), south (i=2), east (i=3), or west (i=4) of the nth island
90      !!      miisl(ip,i,n), mjisl(ip,i,n) : (i,j) index of ipth u- or v-point
91      !!      along (i=0), north (i=1),south (i=2), east (i=3), or west (i=4)
92      !!      of the nth island
93      !!
94      !! History :
95      !!    4.0  !  88-03  (G. Madec)  Original code
96      !!    6.0  !  96-01  (G. Madec)  Suppress common workspace, use of bmask
97      !!    8.5  !  96-01  (G. Madec)  Free form, F90
98      !!----------------------------------------------------------------------
99      !! * Local declarations
100      INTEGER ::   ji, jj, jn, jnil   ! dummy loop indices
101      INTEGER ::   ind, inilt, ip, ipn, ips, ipe, ipw, ii, ij, iju
102      INTEGER ::   iista, iiend, ijsta, ijend, ijstm1, ijenm1
103      INTEGER ::   isrchne
104      INTEGER, DIMENSION(jpj) ::   indil
105      INTEGER, DIMENSION(jpi,jpj) ::   idil
106
107      REAL(wp) znil
108      REAL(wp), DIMENSION(jpi,jpj) ::   zwb
109      !!----------------------------------------------------------------------
110
111
112      ! 0. Islands index computed from mbathy
113      ! -------------------------------------
114      ! zwb=0 over the continent and ocean, =-n over the nth island
115     
116      ! Computation
117      DO jj = 1, jpjm1
118         DO ji = 1, jpim1
119            zwb(ji,jj) = MIN( 0 , mbathy(ji,jj  ), mbathy(ji+1,jj  ),   &
120                                  mbathy(ji,jj+1), mbathy(ji+1,jj+1)  )
121         END DO
122      END DO
123     
124      ! Lateral boundary conditions on zwb
125
126      !  mono- or macro-tasking environnement:
127      IF( nperio == 2 ) THEN
128         zwb(:, 1 ) = zwb(:, 2 )
129      ELSEIF( nperio == 3 .OR. nperio == 4) THEN
130! om
131!$$$          DO ji = 1, jpim1
132! On ne peut pas partir de ji=1, car alors jiu=jpi, et cette valeur
133! n'est pas encore initialiee.
134! Le cas ji=1 est de toute facon traite a partir de la ligne 135
135         DO ji = 2, jpim1
136            iju = jpi-ji+1
137            zwb(ji,jpj  ) = zwb(iju,jpj-3)
138            zwb(ji,jpjm1) = zwb(iju,jpj-2)
139         END DO
140      ELSEIF( nperio == 5 .OR. nperio == 6) THEN
141         DO ji = 1, jpim1
142            iju = jpi-ji
143            zwb(ji,jpj  ) = zwb(iju,jpj-2)
144         END DO
145         DO ji = jpi/2+1, jpi-1
146            iju = jpi-ji
147            zwb(ji,jpjm1) = zwb(iju,jpjm1)
148         END DO
149      ELSE
150         zwb(:,jpj) = 0.e0
151      ENDIF
152
153      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
154         zwb( 1 ,:) = zwb(jpim1,:)
155         zwb(jpi,:) = zwb(  2  ,:)
156      ELSE
157         zwb( 1 ,:) = 0.e0
158         zwb(jpi,:) = 0.e0
159      ENDIF
160      IF( lk_mpp )   CALL lbc_lnk( zwb, 'G', 1. )
161
162
163      ! 1. Initialization for the search of island grid-points
164      ! ------------------------------------------------------
165
166      IF( lk_mpp ) THEN
167
168         ! Mpp : The overlap region are not taken into account
169         ! (islands bondaries are searched over subdomain only)
170         iista =  1   + jpreci
171         iiend = nlci - jpreci
172         ijsta =  1   + jprecj
173         ijend = nlcj - jprecj
174         ijstm1=  1   + jprecj
175         ijenm1= nlcj - jprecj
176         IF( nbondi == -1 .OR. nbondi == 2 ) THEN
177            iista  = 1
178         ENDIF
179         IF( nbondi ==  1 .OR. nbondi == 2 ) THEN
180            iiend  = nlci
181         ENDIF
182         IF( nbondj == -1 .OR. nbondj == 2 ) THEN
183            ijsta  = 1
184            ijstm1 = 2
185         ENDIF
186         IF( nbondj ==  1 .OR. nbondj == 2 ) THEN
187            ijend  = nlcj
188            ijenm1 = nlcj-1
189         ENDIF
190         IF( npolj == 3 .OR. npolj == 4 ) THEN
191            ijend  = nlcj-2
192            ijenm1 = nlcj-2
193         ENDIF
194      ELSE
195         ! mono- or macro-tasking environnement: full domain scan
196         iista  = 1
197         iiend  = jpi
198         ijsta  = 1
199         ijstm1 = 2
200         IF( nperio == 3 .OR. nperio == 4 ) THEN
201            ijend  = jpj-2
202            ijenm1 = jpj-2
203         ELSEIF( nperio == 5 .OR. nperio == 6 ) THEN
204            ijend  = jpj-1
205            ijenm1 = jpj-1
206         ELSE
207            ijend  = jpj
208            ijenm1 = jpj-1
209         ENDIF
210      ENDIF
211
212
213      ! 2. Loop over island
214      ! -------------------
215
216      DO jnil = 1, jpisl
217
218         ! 2.1 Initialization to zero of miisl, mjisl of the jnil island
219
220         miisl(:,:,jnil) = 0
221         mjisl(:,:,jnil) = 0
222         
223         ! 2.2 Search grid-points of island jnil
224 
225         indil(:)   = 0
226         idil (:,:) = 0
227         
228         znil = - FLOAT( jnil )
229         DO jj = ijsta, ijend
230            indil(jj) = 0
231            ind = 0
232            DO ji = iista, iiend
233               IF( zwb(ji,jj) == znil ) THEN
234                  ind = ind + 1
235                  idil(ind,jj) = ji
236                  indil(jj) = indil(jj) + 1
237               ENDIF
238            END DO
239         END DO
240         
241         ! 2.3 Check the number of island
242         
243         inilt = 0
244         DO jj = ijsta, ijend
245            inilt = inilt + indil(jj)
246         END DO
247         IF( lk_mpp )   CALL mpp_sum( inilt )   ! sum over the global domain
248
249         IF( inilt == 0 ) THEN
250            WRITE(ctmp1,*) ' isldom: there is not island number: ', jnil,' while jpisl= ', jpisl
251            CALL ctl_stop( ctmp1, ' change par_oce' )
252
253         ENDIF
254         
255         ! 2.4 Coastal island grid-points (miisl,mjisl)
256         
257         ip  = 0
258         ipn = 0
259         ips = 0
260         ipe = 0
261         ipw = 0
262
263         ! South line (ij=1)
264         ij = 1
265         DO jn = 1, indil(ij)
266            ii = idil(jn,ij)
267            IF( (ij+njmpp-1) == 1 .AND. ii > jpreci .AND. ii < (nlci-jpreci+1) ) THEN
268               IF(  zwb(ii-1,ij) * zwb(ii+1,ij) * zwb(ii,ij+1)  == 0. ) THEN
269                  ip = ip + 1
270                  miisl(ip,0,jnil) = ii
271                  mjisl(ip,0,jnil) = ij
272                  IF( zwb(ii-1,ij) == 0. ) THEN
273                     ipw = ipw + 1
274                     miisl(ipw,4,jnil) = ii
275                     mjisl(ipw,4,jnil) = ij
276                  ENDIF
277                  IF( zwb(ii+1,ij) == 0. ) THEN
278                     ipe = ipe+1
279                     IF( (nperio == 1 .OR. nperio == 4.OR. nperio == 6)  &
280                        .AND. ii == (nlci-jpreci) ) THEN
281                        miisl(ipe,3,jnil) = 1+jpreci
282                     ELSE
283                        miisl(ipe,3,jnil) = ii + 1
284                     ENDIF
285                     mjisl(ipe,3,jnil) = ij
286                  ENDIF
287                  IF( zwb(ii,ij+1) == 0. ) THEN
288                     ipn = ipn+1
289                     miisl(ipn,1,jnil) = ii
290                     mjisl(ipn,1,jnil) = ij + 1
291                  ENDIF
292               ENDIF
293            ENDIF
294         END DO
295         
296         ! Middle lines (2=<jj=<jpjm1 or jpj-2 if north fold b.c.)
297         DO jj = ijstm1, ijenm1
298            DO jn = 1, indil(jj)
299               ii = idil(jn,jj)
300               IF( ii > jpreci .AND. ii < (nlci-jpreci+1) ) THEN
301                  IF( (zwb(ii-1, jj )*zwb(ii+1, jj )*   &
302                       zwb( ii ,jj-1)*zwb( ii ,jj+1)) == 0. ) THEN
303                     ip = ip + 1
304                     miisl(ip,0,jnil) = ii
305                     mjisl(ip,0,jnil) = jj
306                     IF( zwb(ii-1,jj) == 0. ) THEN
307                        ipw = ipw + 1
308                        miisl(ipw,4,jnil) = ii
309                        mjisl(ipw,4,jnil) = jj
310                     ENDIF
311                     IF( zwb(ii+1,jj) == 0. ) THEN
312                        ipe = ipe + 1
313                        IF((nperio == 1.OR.nperio == 4.OR.nperio == 6)   &
314                           .AND.ii == (nlci-jpreci) ) THEN
315                           miisl(ipe,3,jnil) = 1 + jpreci
316                        ELSE
317                           miisl(ipe,3,jnil) = ii + 1
318                        ENDIF
319                        mjisl(ipe,3,jnil) = jj
320                     ENDIF
321                     IF( zwb(ii,jj-1) == 0. ) THEN
322                        ips = ips + 1
323                        miisl(ips,2,jnil) = ii
324                        mjisl(ips,2,jnil) = jj
325                     ENDIF
326                     IF( zwb(ii,jj+1) == 0. ) THEN
327                        ipn = ipn + 1
328                        miisl(ipn,1,jnil) = ii
329                        mjisl(ipn,1,jnil) = jj + 1
330                     ENDIF
331                  ENDIF
332               ENDIF
333            END DO
334         END DO
335         
336         ! North line (jj=jpj) only if not north fold b.c.
337         IF( nperio /= 3 .AND. nperio /= 4 .AND. nperio /= 5 .AND. nperio /= 6 ) THEN
338            ij = jpj
339            DO jn = 1, indil(ij)
340               ii = idil(jn,ij)
341               IF( (ij+njmpp-1) == jpjglo .AND. ii > jpreci .AND. ii < (nlci-jpreci+1) ) THEN
342                  IF( (zwb(ii-1, ij )*zwb(ii+1, ij )* zwb( ii ,ij-1) ) == 0. ) THEN
343                     ip = ip+1
344                     miisl(ip,0,jnil) = ii
345                     mjisl(ip,0,jnil) = ij
346                     IF( zwb(ii-1,ij) == 0. ) THEN
347                        ipw = ipw+1
348                        miisl(ipw,4,jnil) = ii
349                        mjisl(ipw,4,jnil) = ij
350                     ENDIF
351                     IF( zwb(ii+1,ij) == 0. ) THEN
352                        ipe = ipe+1
353                        IF( (nperio == 1) .AND. ii == (nlci-jpreci) ) THEN
354                           miisl(ipe,3,jnil) = 1+jpreci
355                        ELSE
356                           miisl(ipe,3,jnil) = ii+1
357                        ENDIF
358                        mjisl(ipe,3,jnil) = ij
359                     ENDIF
360                     IF( zwb(ii,ij-1) == 0. ) THEN
361                        ips = ips+1
362                        miisl(ips,2,jnil) = ii
363                        mjisl(ips,2,jnil) = ij
364                     ENDIF
365                  ENDIF
366               ENDIF
367            END DO
368         ENDIF
369         
370         mnisl(0,jnil) = ip
371         mnisl(1,jnil) = ipn
372         mnisl(2,jnil) = ips
373         mnisl(3,jnil) = ipe
374         mnisl(4,jnil) = ipw
375         
376         ! Take account of redundant points
377         
378         IF( lk_mpp )   CALL mpp_sum( ip )   ! sum over the global domain
379         
380         IF( ip > jpnisl ) THEN
381            WRITE(ctmp1,*) ' isldom: the island ',jnil,' has ',   &
382                 mnisl(0,jnil),' grid-points, while jpnisl= ', jpnisl,ip
383            CALL ctl_stop( ctmp1, ' change par_oce.h' )
384         ENDIF
385         
386         ! 2.5 Set to zero the grid-points of the jnil island in x
387         
388         DO jj = 1, jpj
389            DO ji = 1, jpi
390               IF( zwb(ji,jj)+FLOAT(jnil) == 0. ) zwb(ji,jj) = 0. 
391            END DO
392         END DO
393         
394      END DO
395     
396
397      ! 3. Check the number of island
398      ! -----------------------------
399
400      inilt = isrchne( jpij, zwb(1,1), 1, 0. )
401      IF( lk_mpp )   CALL mpp_min( inilt )   ! min over the global domain
402
403      IF( inilt /= jpij+1 ) THEN
404            WRITE(ctmp1,*) ' isldom: there is at least one more ',   &
405                  'island in the domain and jpisl=', jpisl
406            CALL ctl_stop( ctmp1, ' change par_oce.h' )
407      ENDIF
408
409
410      ! 4. Print of island parametres and arrays
411      ! ----------------------------------------
412
413      CALL isl_pri
414
415
416      ! 5. Array for computation of circulation arround islands
417      ! -------------------------------------------------------
418
419      CALL isl_pth
420
421   END SUBROUTINE isl_dom
422
423
424   SUBROUTINE isl_pri
425      !!----------------------------------------------------------------------
426      !!                  ***  ROUTINE isl_pri  ***
427      !!             
428      !! ** Purpose :   Print islands variables and islands arrays
429      !!
430      !! ** Method  :
431      !!
432      !! History :
433      !!   1.0  !  88-03  (G. Madec)  Original code
434      !!   8.5  !  02-08  (G. Madec)  Free form, F90
435      !!----------------------------------------------------------------------
436      !! * Local declarations
437      INTEGER ::   ji, jni, jnp      ! dummy loop variables
438      INTEGER ::   &
439         ip, ipn, ips, ipe, ipw      ! temporary integers
440      !!----------------------------------------------------------------------
441     
442      IF(lwp) THEN
443         WRITE(numout,*) ' '
444         WRITE(numout,*) '*** islpri number of islands : jpisl =',jpisl
445      ENDIF
446
447      DO jni = 1, jpisl
448         ip  = mnisl(0,jni)
449         ipn = mnisl(1,jni)
450         ips = mnisl(2,jni)
451         ipe = mnisl(3,jni)
452         ipw = mnisl(4,jni)
453         IF( lk_mpp ) THEN
454            CALL mpp_sum( ip  )   ! sums over the global domain
455            CALL mpp_sum( ipn )
456            CALL mpp_sum( ips )
457            CALL mpp_sum( ipe )
458            CALL mpp_sum( ipw )
459         ENDIF
460         IF(lwp) THEN
461            WRITE(numout,9000) jni
462            WRITE(numout,9010) ip, ipn, ips, ipe, ipw
463            WRITE(numout,9020)
464            DO jnp = 1, mnisl(0,jni)
465               WRITE(numout,9030) jnp, ( miisl(jnp,ji,jni)+nimpp-1,   &
466                                         mjisl(jnp,ji,jni)+njmpp-1, ji = 0, 4 )
467            END DO
468         ENDIF
469      END DO
470
471      ! FORMAT   !!cr => no more format
472 9000 FORMAT(/, /, 'island number= ', i2 )
473 9010 FORMAT(/, 'npil=',i4,' npn=',i3,' nps=',i3,' npe=',i3,' npw=',i3 )
474 9020 FORMAT(/,'     * isl point *  point n  *  point s  *  point e ', '*  point w  *')
475 9030 FORMAT(i4,' * (',i3,',',i3,') * (',i3,',',i3,') * (',i3,',',i3,   &
476          ') * (',i3,',',i3,') * (',i3,',',i3,') *' )
477
478   END SUBROUTINE isl_pri
479
480
481   SUBROUTINE isl_pth
482      !!----------------------------------------------------------------------
483      !!                  ***  ROUTINE isl_pth  ***
484      !!             
485      !! ** Purpose :   intialize arrays for the computation of streamfunction
486      !!      around islands
487      !!
488      !! ** Method  :
489      !!
490      !! ** Action  : - acisl1(i,j,ii,n): coefficient n-s (ii=1) and e-w (ii=2)
491      !!                for the calculation of mu and mv around the island n     
492      !!              - acisl2(i,j,ii,n): coefficient n-s (ii=1) and e-w (ii=2)
493      !!                for the calculation of bsfd around the island n
494      !!
495      !! History :
496      !!   1.0  !  88-03  (G. Madec)  Original code
497      !!   8.5  !  02-08  (G. Madec)  Free form, F90
498      !!----------------------------------------------------------------------
499      !! * Local declarations
500      INTEGER ::   jni, jii, jnp    ! dummy loop indices
501      INTEGER ::   ii, ij           ! temporary integers
502      !!----------------------------------------------------------------------
503     
504      ! 1. Initialisation
505      ! -----------------
506      acisl1(:,:,:,:) = 0.e0
507      acisl2(:,:,:,:) = 0.e0
508      bsfisl(:,:,:)   = 0.e0
509
510      ! 2. Coefficient arrays
511      ! ---------------------
512      DO jni = 1, jpisl
513         DO jii = 1, 4
514            DO jnp = 1,mnisl(jii,jni)
515               ii = miisl(jnp,jii,jni)
516               ij = mjisl(jnp,jii,jni)
517               IF( jii <= 2 ) THEN
518                  ! north and south points
519                  acisl1(ii,ij,1,jni) = float( 2*jii-3) * e1u(ii,ij)
520                  acisl2(ii,ij,1,jni) = float( 2*jii-3) * e1u(ii,ij) * hur(ii,ij) / e2u(ii,ij)
521               ELSE
522                  ! east and west points
523                  acisl1(ii,ij,2,jni) = float(-2*jii+7) * e2v(ii,ij)
524                  acisl2(ii,ij,2,jni) = float(-2*jii+7) * e2v(ii,ij) * hvr(ii,ij) / e1v(ii,ij)
525               ENDIF
526            END DO
527         END DO
528      END DO
529     
530   END SUBROUTINE isl_pth
531
532   !!----------------------------------------------------------------------
533   !!   Default option :                                        NetCDF file
534   !!----------------------------------------------------------------------
535
536   SUBROUTINE isl_mat
537      !!----------------------------------------------------------------------
538      !!                  ***  ROUTINE isl_mat  ***
539      !!               
540      !! ** Purpose :   Compute and invert the island matrix
541      !!
542      !! ** Method  :   aisl(jni,jnj) matrix constituted by bsfisl circulation
543      !!
544      !! ** Action  : - aisl     : island matrix
545      !!              - aislm1   : invert of the island matrix
546      !!      file
547      !!              - islands  : contain bsfisl, aisl and aislm1
548      !!
549      !! History :
550      !!   1.0  !  88-10  (G. Madec)  Original code
551      !!   7.0  !  96-01  (G. Madec)  suppression of common work arrays
552      !!   8.5  !  02-08  (G. Madec)  Free form, F90
553      !!----------------------------------------------------------------------
554      !! * Modules used
555      USE ioipsl
556      USE iom
557     
558      !! * Local declarations
559      INTEGER ::   ji, jj, jni, jnj, jl   ! dummy loop indices
560      INTEGER ::   ios          ! temporary integers
561      INTEGER  ::   &
562         inum                 ! temporary logical unit
563      REAL(wp), DIMENSION(jpi,jpj,2) ::   zwx
564      REAL(wp), DIMENSION(jpisl*jpisl) ::   ztab
565      !!----------------------------------------------------------------------
566
567
568      ! I. Island matrix lecture in numisl (if it exists)
569      ! ==================================
570
571      ! Lecture
572      CALL iom_open ( 'islands', inum )
573      ios = iom_varid( inum, 'aisl', ldstop = .FALSE. )
574      IF( ios > 0 ) THEN
575
576         CALL iom_get( inum, jpdom_unknown, 'aisl'  , aisl )     
577         CALL iom_get( inum, jpdom_unknown, 'aislm1', aislm1 )     
578         CALL iom_close( inum )
579         ! Control print
580         IF(lwp) THEN
581            WRITE(numout,*)
582            WRITE(numout,*)' islmat: lecture aisl/aislm1 in numisl done'
583            WRITE(numout,*)' ~~~~~~'
584            WRITE(numout,*)
585            WRITE(numout,*) '        island matrix : '
586            WRITE(numout,*)
587           
588            DO jnj = 1, jpisl
589               WRITE(numout,'(8e12.4)') ( aisl(jni,jnj), jni = 1, jpisl )
590            END DO
591           
592            WRITE(numout,*)
593            WRITE(numout,*) '       inverse of the island matrix'
594            WRITE(numout,*)
595           
596            DO jnj = 1, jpisl
597               WRITE(numout,'(12e11.3)') ( aislm1(jni,jnj), jni=1,jpisl )
598            END DO
599         ENDIF
600 
601         CALL restclo(numisl)
602                 
603      ELSE
604 
605         CALL iom_close( inum )
606
607      ! II. Island matrix computation
608      ! =============================
609
610      DO jnj = 1, jpisl
611
612         ! Circulation of bsf(jnj) around island jni
613         
614         DO jj = 2, jpj
615               zwx(:,jj,1) = -( bsfisl(:,jj,jnj) - bsfisl(:,jj-1,jnj) )
616         END DO
617         zwx(:,1,1) = 0.e0
618         
619         DO jj = 1, jpj
620            DO ji = 2, jpi
621               zwx(ji,jj,2) = ( bsfisl(ji,jj,jnj) - bsfisl(ji-1,jj,jnj) )
622            END DO
623         END DO
624         zwx(1,:,2) = 0.e0
625         
626         ! Island matrix
627         
628         DO jni = 1, jpisl
629            aisl(jni,jnj) = 0.e0
630            DO jl = 1, 2
631               DO jj=1,jpj
632                  DO ji=1,jpi
633                     aisl(jni,jnj) = aisl(jni,jnj) + acisl2(ji,jj,jl,jni)*zwx(ji,jj,jl)
634                  END DO
635               END DO
636            END DO
637         END DO
638         
639      END DO
640      IF( lk_mpp ) THEN
641         DO jnj = 1, jpisl
642            DO jni = 1, jpisl
643               ztab(jni+(jnj-1)*jpisl) = aisl(jni,jnj)
644            END DO
645         END DO
646         CALL mpp_sum( ztab, jpisl*jpisl )   ! sum over the global domain
647      ENDIF
648
649      ! 1.3 Control print
650
651      IF(lwp) THEN
652         WRITE(numout,*)
653         WRITE(numout,*) 'islmat : island matrix'
654         WRITE(numout,*) '~~~~~~'
655         WRITE(numout,*)
656         
657         DO jnj = 1, jpisl
658            WRITE(numout,'(8e12.4)') ( aisl(jni,jnj), jni = 1, jpisl )
659         END DO
660      ENDIF
661     
662
663      ! 2. Invertion of the island matrix
664      ! ---------------------------------
665     
666      ! 2.1 Call of an imsl routine for the matrix invertion
667
668      CALL linrg( jpisl, aisl, jpisl, aislm1, jpisl )
669
670      ! 2.2 Control print
671     
672      IF(lwp) THEN
673         WRITE(numout,*)
674         WRITE(numout,*) 'islmat : inverse of the island matrix'
675         WRITE(numout,*) '~~~~~~'
676         WRITE(numout,*)
677         
678         DO jnj = 1, jpisl
679            WRITE(numout, '(12e11.3)')  '        ', ( aislm1(jni,jnj), jni=1, jpisl )
680         END DO
681      ENDIF
682     
683
684      ! 3. Output of aisl and aislm1 in numisl
685      ! --------------------------------------
686
687      CALL restput( numisl, 'aisl'  , jpisl, jpisl, 1, 0, aisl   )
688      CALL restput( numisl, 'aislm1', jpisl, jpisl, 1, 0, aislm1 )
689      CALL restclo( numisl )
690
691      ENDIF
692
693   END SUBROUTINE isl_mat
694
695
696   SUBROUTINE isl_bsf
697      !!----------------------------------------------------------------------
698      !!                  ***  ROUTINE isl_bsf  ***
699      !!         
700      !! ** Purpose :
701      !!         Compute the barotropic stream function associated with each
702      !!      island using FETI or preconditioned conjugate gradient method
703      !!      or read them in numisl.
704      !!
705      !! ** Method :
706      !!
707      !! ** input/output file :
708      !!            numisl            : barotropic stream function associated
709      !!                          with each island of the domain
710      !!
711      !! ** Action :
712      !!       bsfisl, the streamfunction which takes the value 1 over island ni,
713      !!    and 0 over the others islands
714      !!       file 'numisl' barotropic stream function associated
715      !!                              with each island of the domain
716      !!
717      !! History :
718      !!        !  87-10  (G. Madec)  Original code
719      !!        !  91-11  (G. Madec)
720      !!        !  93-03  (G. Madec)  release 7.1
721      !!        !  93-04  (M. Guyon)  loops and suppress pointers
722      !!        !  96-11  (A. Weaver)  correction to preconditioning
723      !!        !  98-02  (M. Guyon)  FETI method
724      !!        !  99-11  (M. Imbard)  NetCDF FORMAT with IOIPSL
725      !!   8.5  !  02-08  (G. Madec)  Free form, F90
726      !!----------------------------------------------------------------------
727      !! * Modules used
728      USE ioipsl
729      USE iom
730      USE solpcg
731      USE solfet
732      USE solsor
733
734      !! * Local declarations
735      LOGICAL  ::   llog, llbon
736      CHARACTER (len=10) ::  clisl
737      CHARACTER (len=32) ::  clname = 'islands'
738      INTEGER  ::   &
739         inum                 ! temporary logical unit
740      INTEGER  ::   ji, jj, jni, jii, jnp  ! dummy loop indices
741      INTEGER  ::   iimlu, ijmlu, inmlu, iju
742      INTEGER  ::   ii, ij, icile, icut, inmax, indic
743      INTEGER  ::   itime
744      REAL(wp) ::   zepsr, zeplu, zgwgt
745      REAL(wp) ::   zep(jpisl), zdept(1), zprec(4)
746      REAL(wp) ::   zdate0, zdt
747      REAL(wp) ::   t2p1(jpi,1,1)
748      INTEGER  ::   iloc
749      !!----------------------------------------------------------------------
750
751
752      ! 0. Initializations
753      ! ==================
754
755      icile         = 0      ! set to zero the convergence indicator
756      bsfisl(:,:,:) = 0.e0   ! set to zero of bsfisl
757
758
759      ! I. Lecture of bsfisl in numisl (if it exists)
760      ! =============================================
761
762      icut  = 0
763      iimlu = 0
764      ijmlu = 0
765      inmlu = 0
766      zeplu = 0.
767     
768      clname = 'islands'
769     
770      INQUIRE( FILE=clname, EXIST=llbon )
771! islands FILE does not EXIST : icut=999
772      IF( llbon ) THEN 
773
774         ! island FILE is present
775
776         CALL iom_open (clname, inum )
777         CALL iom_get( inum, jpdom_unknown, 'PRECISION', zprec )     
778
779         iimlu = NINT( zprec(1) )
780         ijmlu = NINT( zprec(2) )
781         inmlu = NINT( zprec(3) )
782         zeplu = zprec(4)
783         ! the read domain does not correspond to the model one : icut=999
784         IF( iimlu /= jpi .OR. ijmlu /= jpj .OR. inmlu /= jpisl ) THEN
785            icut = 999
786         ELSE
787            DO jni = 1, jpisl
788               IF( jni < 10 ) THEN
789                  WRITE(clisl,'("island",I1)') jni
790               ELSEIF( jni < 100 ) THEN
791                  WRITE(clisl,'("island",I2)') jni
792               ELSE
793                  WRITE(clisl,'("island",I3)') jni
794               ENDIF
795               CALL iom_get( inum, jpdom_autoglo, clisl, bsfisl(:,:,jni))     
796            END DO
797         ENDIF
798      ELSE 
799         ! islands FILE does not EXIST : icut=999
800         icut = 999
801      ENDIF
802 
803      CALL iom_close( inum )
804   
805      ! the read precision is not the required one : icut=888
806      IF( zeplu > epsisl ) THEN
807         icut = 888
808      ENDIF
809
810      ! Control print
811      IF( icut == 999 ) THEN
812          IF(lwp) THEN
813              WRITE(numout,*)
814              WRITE(numout,*) 'islbsf : lecture bsfisl in numisl failed'
815              WRITE(numout,*) '~~~~~~'
816              WRITE(numout,*) '         icut= ', icut
817              WRITE(numout,*) '         imlu= ', iimlu, ' jmlu= ', ijmlu,   &
818                                      ' nilu= ', inmlu, ' epsisl lu= ', zeplu
819              WRITE(numout,*) '         the bsfisl are computed from zero'
820          ENDIF
821        ELSEIF( icut == 888 ) THEN
822          IF(lwp) THEN
823              WRITE(numout,*)
824              WRITE(numout,*) 'islbsf : lecture bsfisl in numisl done'
825              WRITE(numout,*) '~~~~~~'
826              WRITE(numout,*) '         the required accuracy is not reached'
827              WRITE(numout,*) '         epsisl lu= ', zeplu,' epsisl required ', epsisl
828              WRITE(numout,*) '         the bsfisl are computed from the read values'
829          ENDIF
830        ELSE
831          IF(lwp) THEN
832              WRITE(numout,*)
833              WRITE(numout,*) 'islbsf : lecture bsfisl in numisl done'
834              WRITE(numout,*) '~~~~~~'
835          ENDIF
836          RETURN
837      ENDIF
838     
839     
840      ! II. Compute the bsfisl (if icut=888 or 999)
841      ! ============================================
842
843      ! save nmax
844      inmax = nmax
845
846      ! set the number of iteration of island computation
847      nmax = nmisl
848
849      ! Loop over islands
850      ! -----------------
851     
852      DO jni = 1, jpisl
853
854
855         ! 1. Initalizations of island computation
856         ! ---------------------------------------
857         
858         ! Set the pcg solution gcb to zero
859         gcb(:,:) = 0.e0
860
861         ! Set first guess gcx either to zero or to the read bsfisl
862         IF( icut == 999 ) THEN
863            gcx(:,:) = 0.e0
864         ELSEIF( icut == 888 ) THEN
865            IF(lwp) WRITE(numout,*) ' islbsf: bsfisl read are used as first guess'
866            ! c a u t i o n: bsfisl masked because it contains 1 along island seaside
867            gcx(:,:) = bsfisl(:,:,jni) * bmask(:,:)
868         ENDIF
869         
870         ! Right hand side of the streamfunction equation
871         
872         IF( lk_mpp ) THEN
873
874            ! north fold treatment
875            IF( npolj == 3 .OR. npolj == 5)   iloc=jpiglo-(nimpp-1+nimppt(nono+1)-1)
876            IF( npolj == 4 .OR. npolj == 6)   iloc=jpiglo-2*(nimpp-1)
877            t2p1(:,1,1) = 0.e0
878            ! north and south grid-points
879            DO jii = 1, 2
880               DO jnp = 1, mnisl(jii,jni)
881                  ii = miisl(jnp,jii,jni)
882                  ij = mjisl(jnp,jii,jni)
883                  IF( ( npolj == 3 .OR. npolj == 4 ) .AND.   &
884                     ( ij == nlcj-1 .AND. jii == 1) ) THEN
885                     iju=iloc-ii+1
886                     t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
887                  ELSEIF( ( npolj == 5 .OR. npolj == 6 ) .AND.   &
888                     ( ij == nlcj-1 .AND. jii == 1) ) THEN
889                     iju=iloc-ii
890                     gcb(ii,ij) =  gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
891                     t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
892                  ELSE 
893                     gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
894                  ENDIF
895               END DO
896            END DO
897         
898            ! east and west grid-points
899         
900            DO jii = 3, 4
901               DO jnp = 1, mnisl(jii,jni)
902                  ii = miisl(jnp,jii,jni)
903                  ij = mjisl(jnp,jii,jni)
904                  gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
905               END DO
906            END DO
907
908            IF( lk_mpp )   CALL mpplnks( gcb )   !!bug ? should use an lbclnk ? is it possible?
909
910         ELSE
911
912            ! north and south grid-points
913            DO jii = 1, 2
914               DO jnp = 1, mnisl(jii,jni)
915                  ii = miisl(jnp,jii,jni)
916                  ij = mjisl(jnp,jii,jni)
917                  IF( ( nperio == 3 .OR. nperio == 4 ) .AND.   &
918                     ( ij == jpj-1 .AND. jii == 1) ) THEN
919                     gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
920                  ELSEIF( ( nperio == 5 .OR. nperio == 6 ) .AND.   &
921                     ( ij == jpj-1 .AND. jii == 1) ) THEN
922                     gcb(ii,ij) =  gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
923                     gcb(jpi-ii,ij) = gcb(jpi-ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
924                  ELSE 
925                     gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
926                  ENDIF
927               END DO
928            END DO
929
930            ! east and west grid-points
931            DO jii = 3, 4
932               DO jnp = 1, mnisl(jii,jni)
933                  ii = miisl(jnp,jii,jni)
934                  ij = mjisl(jnp,jii,jni)
935                  IF( bmask(ii-jii+3,ij) /= 0. ) THEN
936                     gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
937                  ELSE
938                     ! east-west cyclic boundary conditions
939                     IF( ii-jii+3 == 1 ) THEN
940                        gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
941                     ENDIF
942                  ENDIF
943               END DO
944            END DO
945         ENDIF
946
947         ! Preconditioned right hand side and absolute precision
948
949         IF( nsolv == 3 ) THEN 
950            ! FETI method
951            ncut    = 0
952            rnorme  = 0.e0
953            gcb(:,:) = bmask(:,:) * gcb(:,:)
954            DO jj = 1, jpj
955               DO ji = 1, jpi
956                  rnorme = rnorme + gcb(ji,jj) * gcb(ji,jj)
957               END DO
958            END DO
959           
960            IF( lk_mpp )   CALL mpp_sum( rnorme )
961
962            IF(lwp) WRITE(numout,*) 'rnorme ', rnorme
963            epsr  = epsisl * epsisl * rnorme
964            indic = 0
965         ELSE
966            ncut   = 0
967            rnorme = 0.e0
968            DO jj = 1, jpj
969               DO ji = 1, jpi
970                  gcb  (ji,jj) = gcdprc(ji,jj) * gcb(ji,jj)
971                  zgwgt  = gcdmat(ji,jj) * gcb(ji,jj)
972                  rnorme = rnorme + gcb(ji,jj) * zgwgt
973               END DO
974            END DO
975            IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
976
977            IF(lwp) WRITE(numout,*) 'rnorme ', rnorme
978            epsr = epsisl * epsisl * rnorme
979            indic = 0
980         ENDIF
981
982
983         ! 3. PCG solver for gcp.gcx=gcb in monotask
984         ! -------------------------------------------
985
986         IF( nsolv == 3 ) THEN
987            epsilo = epsisl       ! precision to compute Islands matrix A
988            CALL sol_fet( indic )  ! FETI method
989            epsilo = eps          ! precision to compute grad PS
990         ELSE
991            CALL sol_pcg( indic )  ! pcg method
992         ENDIF
993
994
995         ! 4. Save the solution in bsfisl
996         ! ------------------------------
997
998         bsfisl(:,:,jni) = gcx(:,:)
999
1000
1001         ! 5. Boundary conditions
1002         ! ----------------------
1003
1004         ! set to 1. coastal gridpoints of the island
1005         DO jnp = 1, mnisl(0,jni)
1006            ii = miisl(jnp,0,jni)
1007            ij = mjisl(jnp,0,jni)
1008            bsfisl(ii,ij,jni) = 1.
1009         END DO
1010
1011         ! cyclic boundary conditions
1012         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
1013            bsfisl( 1 ,:,jni) = bsfisl(jpim1,:,jni)
1014            bsfisl(jpi,:,jni) = bsfisl(  2  ,:,jni)
1015         ENDIF
1016         IF( nperio == 3 .OR. nperio == 4 ) THEN
1017            DO ji = 1, jpim1
1018               iju = jpi-ji+1
1019               bsfisl(ji,jpj-1,jni) = bsfisl(iju,jpj-2,jni)
1020               bsfisl(ji, jpj ,jni) = bsfisl(iju,jpj-3,jni)
1021            END DO
1022         ENDIF
1023         IF( nperio == 5 .OR. nperio == 6 ) THEN
1024            DO ji = 1, jpi-1
1025               iju=jpi-ji
1026               bsfisl(ji,jpj,jni) = bsfisl(iju,jpj-2,jni)
1027            END DO
1028            DO ji = jpi/2+1, jpi-1
1029               iju=jpi-ji
1030               bsfisl (ji,jpjm1,jni) = bsfisl(iju,jpjm1,jni)
1031            END DO
1032         ENDIF
1033         IF( lk_mpp )   CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. )   ! link at G-point
1034
1035
1036         ! 6. Control print
1037         ! ----------------
1038
1039         IF(lwp) WRITE(numout,*)
1040         IF(lwp) WRITE(numout,*) ' islbsf: island number: ', jni
1041         IF(lwp) WRITE (numout,9290) niter, res, SQRT(epsr)/epsisl
1042         zep(jni) =  MAX(epsisl, res/(SQRT(epsr)/epsisl))
1043         IF( indic <  0 ) THEN
1044            icile = icile-1
1045            IF(lwp) WRITE(numout,*) '  pcg do not converge for island: ', jni
1046            IF(lwp) WRITE(numout,*) '      Precision reached: ',zep(jni)
1047         ENDIF
1048
10499290     FORMAT('        niter :',i4,' , res :',e20.10,' , gcb :',e20.10)
1050
1051         !                                          !====================
1052      END DO                                        !  End Loop islands
1053      !                                             !====================
1054
1055      ! 7. Reset PCG
1056      ! ------------
1057
1058      ! reset the number of iteration for pcg
1059      nmax = inmax
1060
1061      ! reset to zero pcg arrays
1062      gcx  (:,:) = 0.e0
1063      gcxb (:,:) = 0.e0
1064      gcb  (:,:) = 0.e0
1065      gcr  (:,:) = 0.e0
1066      gcdes(:,:) = 0.e0
1067      gccd (:,:) = 0.e0
1068
1069
1070      ! III. Output of bsfisl in numisl
1071      ! ===============================
1072
1073      CALL ymds2ju( 0, 1, 1, 0.e0, zdate0 )
1074      zprec(1) = FLOAT(jpi)
1075      zprec(2) = FLOAT(jpj)
1076      zprec(3) = FLOAT(jpisl)
1077      IF(lwp) WRITE(numout,*) clname
1078      zdept(1)   = 0.
1079      itime = 0
1080      CALL restini( 'NONE', jpi, jpj, glamt, gphit, 1, zdept,  &
1081         &          clname, itime, zdate0, rdt, numisl, domain_id=nidom )
1082      IF( icile == 0 .AND. icut /= 0 ) THEN
1083         IF(lwp) THEN
1084            WRITE(numout,*)
1085            WRITE(numout,*)' islbsf: write bsfisl in numisl ', numisl
1086            WRITE(numout,*)' -------------------------------------'
1087         ENDIF
1088         zprec(4) = epsisl
1089         CALL restput(numisl,'PRECISION',1,1,4,0,zprec)
1090         DO jni = 1, jpisl
1091            IF(jni < 10) THEN
1092               WRITE(clisl,'("island",I1)') jni
1093            ELSE IF(jni < 100) THEN
1094               WRITE(clisl,'("island",I2)') jni
1095            ELSE
1096               WRITE(clisl,'("island",I3)') jni
1097            ENDIF
1098            CALL restput( numisl, clisl, jpi, jpj, 1, 0, bsfisl(:,:,jni) )
1099         END DO
1100      ENDIF
1101
1102      IF( icile < 0 ) THEN
1103         IF(lwp) THEN
1104            WRITE(numout,*)
1105            WRITE(numout,*) ' islbsf: number of island without convergence : ',ABS(icile)
1106            WRITE(numout,*) ' ---------------------------------------------'
1107         ENDIF
1108         zepsr = epsisl
1109         DO jni = 1, jpisl
1110            IF(lwp) WRITE(numout,*) '    isl ',jni,' precision reached ', zep(jni)
1111            zepsr = MAX( zep(jni), zepsr )
1112         END DO
1113         IF( zepsr == 0. ) zepsr = epsisl
1114         IF(lwp) THEN
1115            WRITE(numout,*) ' save value of precision reached: ',zepsr
1116            WRITE(numout,*)
1117            WRITE(numout,*)' islbsf: save bsfisl in numisl ',numisl
1118            WRITE(numout,*)' -------------------------------------'
1119         ENDIF
1120
1121         zprec(4) = zepsr
1122         CALL restput( numisl, 'PRECISION', 1, 1, 1, 0, zprec )
1123         DO jni = 1, jpisl
1124            IF( jni < 10 ) THEN
1125               WRITE(clisl,'("island",I1)') jni
1126            ELSE IF( jni < 100 ) THEN
1127               WRITE(clisl,'("island",I2)') jni
1128            ELSE
1129               WRITE(clisl,'("island",I3)') jni
1130            ENDIF
1131            CALL restput( numisl, clisl, jpi, jpj, 1, 0, bsfisl(:,:,jni) )
1132         END DO
1133         CALL restclo(numisl)
1134         CALL ctl_stop( ' ' )
1135      ENDIF
1136
1137   END SUBROUTINE isl_bsf
1138
1139
1140   SUBROUTINE isl_dyn_spg
1141      !!----------------------------------------------------------------------
1142      !!                  ***  routine isl_dyn_spg  ***
1143      !!
1144      !! ** Purpose :   Compute and add the island contribution to the
1145      !!      barotropic stream function trend.
1146      !!
1147      !!
1148      !! ** Method  :   Rigid-lid appromimation: ...????
1149      !!
1150      !! ** Action : - Update bsfd with the island contribution
1151      !!
1152      !! History :
1153      !!   9.0  !  03-09  (G. Madec)  isolate island computation
1154      !!---------------------------------------------------------------------
1155
1156      !! * Local declarations
1157      INTEGER ::   ji, jj, jni, jnj    ! dummy loop indices
1158      !!----------------------------------------------------------------------
1159     
1160
1161      ! compute the island potential
1162      ! ----------------------------
1163      DO jni = 1, jpisl                      ! second member
1164         bisl(jni) =  0.e0
1165         DO jj = 2, jpj
1166            DO ji = 2, jpi
1167               bisl(jni) =  bisl(jni) + acisl1(ji,jj,1,jni) *  spgu(ji,jj)                  &
1168                  &                   + acisl1(ji,jj,2,jni) *  spgv(ji,jj)                  &
1169                  &                   + acisl2(ji,jj,1,jni) * ( gcx(ji,jj)-gcx(ji,jj-1) )   &
1170                  &                   - acisl2(ji,jj,2,jni) * ( gcx(ji,jj)-gcx(ji-1,jj) )
1171            END DO
1172         END DO
1173      END DO
1174      IF( lk_mpp )   CALL mpp_sum( bisl, jpisl )   ! sum over the global domain
1175
1176      DO jni = 1, jpisl                     ! Island stream function trend
1177         visl(jni) = 0.e0
1178         DO jnj = 1, jpisl
1179            visl(jni) = visl(jni) + aislm1(jni,jnj) * bisl(jnj)
1180         END DO
1181      END DO
1182
1183      ! update the bsf trend ( caution : bsfd is not zero along island coastlines, dont mask it ! )
1184      ! --------------------
1185      DO jj = 1, jpj
1186         DO jni = 1, jpisl
1187            DO ji = 1, jpi
1188               bsfd(ji,jj) = bsfd(ji,jj) + visl(jni) * bsfisl(ji,jj,jni)
1189            END DO
1190         END DO
1191      END DO
1192
1193   END SUBROUTINE isl_dyn_spg
1194
1195
1196   SUBROUTINE isl_stp_ctl( kt, kindic )
1197      !!----------------------------------------------------------------------
1198      !!                    ***  ROUTINE isl_stp_ctl  ***
1199      !!                     
1200      !! ** Purpose :   ???
1201      !!
1202      !! ** Method  : - print island potential
1203      !!
1204      !! History :
1205      !!   9.0  !  03-09  (G. Madec)  isolated from stp_ctl
1206      !!----------------------------------------------------------------------
1207      !! * Arguments
1208      INTEGER, INTENT(in   ) ::   kt        ! ocean time-step index
1209      INTEGER, INTENT(inout) ::   kindic    ! indicator of solver convergence
1210
1211      !! * local declarations
1212      INTEGER  ::   jni                     ! dummy loop indice
1213      REAL(wp) ::   zfact                   ! temporary scalar
1214      !!----------------------------------------------------------------------
1215      !!  OPA 8.5, LODYC-IPSL (2002)
1216      !!----------------------------------------------------------------------
1217
1218      IF( kt == nit000 .AND. lwp ) THEN
1219         WRITE(numout,*)
1220         WRITE(numout,*) 'isl_stp_ctl : time-stepping control'
1221         WRITE(numout,*) '~~~~~~~~~~~'
1222      ENDIF
1223
1224      ! Island trends
1225      DO jni = 1, jpisl
1226         zfact = 0.
1227         IF( miisl(1,0,jni) /= 0 .AND. mjisl(1,0,jni) /= 0 ) THEN
1228            zfact = 1.e-6 * bsfn(miisl(1,0,jni),mjisl(1,0,jni))
1229         ENDIF
1230         IF( lk_mpp )   CALL mpp_isl( zfact )
1231
1232         IF(lwp) WRITE(numisp,9300) kt, jni, zfact, visl(jni)
1233         IF( MOD( kt, nwrite ) == 0 .OR. kindic < 0     &
1234            .OR. ( kt == nit000 .AND. kindic > 0 )      &
1235            .OR.   kt == nitend                         ) THEN
1236            IF( jni == 1 .AND. lwp ) THEN
1237               WRITE(numout,*)
1238               WRITE(numout,*) 'isl_stp_ctl : island bsf'
1239               WRITE(numout,*) '~~~~~~~~~~~'
1240            ENDIF
1241            IF(lwp) WRITE(numout,9300) kt, jni, zfact, visl(jni)
1242         ENDIF
1243      END DO
12449300  FORMAT(' it : ',i8,' island  :',i4,'   BSF (Sverdrup) : ',f7.2, '   visl : ',e15.6)
1245
1246   END SUBROUTINE isl_stp_ctl
1247
1248#else
1249   !!----------------------------------------------------------------------
1250   !!   Default option                                         Empty module
1251   !!----------------------------------------------------------------------
1252   LOGICAL, PUBLIC, PARAMETER ::   lk_isl = .FALSE.    !: 'key_islands' flag
1253CONTAINS
1254   SUBROUTINE isl_dom                        ! Empty routine
1255   END SUBROUTINE isl_dom
1256   SUBROUTINE isl_bsf                        ! Empty routine
1257   END SUBROUTINE isl_bsf
1258   SUBROUTINE isl_mat                        ! Empty routine
1259   END SUBROUTINE isl_mat
1260   SUBROUTINE isl_dyn_spg                    ! Empty routine
1261   END SUBROUTINE isl_dyn_spg
1262   SUBROUTINE isl_stp_ctl( kt, kindic )      ! Empty routine
1263      WRITE(*,*) 'isl_stp_ctl: You should not have seen this print! error?', kt, kindic
1264   END SUBROUTINE isl_stp_ctl
1265#endif
1266
1267   !!======================================================================
1268END MODULE solisl
Note: See TracBrowser for help on using the repository browser.