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.
solisl.F90 in trunk/NEMO/OPA_SRC/SOL – NEMO

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

Last change on this file since 1473 was 1415, checked in by ctlod, 15 years ago

replace obsolete restini, restput & restclo using IOM subroutines, see ticket: #423

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 44.3 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   !! $Id$
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 ::   inum                   ! temporary logical unit
562      REAL(wp), DIMENSION(jpi,jpj,2) ::   zwx
563      REAL(wp), DIMENSION(jpisl*jpisl) ::   ztab
564      !!----------------------------------------------------------------------
565
566
567      ! I. Island matrix lecture in numisl (if it exists)
568      ! ==================================
569
570      ! Lecture
571      CALL iom_open ( 'islands', inum )
572      ios = iom_varid( inum, 'aisl', ldstop = .FALSE. )
573      IF( ios > 0 ) THEN
574
575         CALL iom_get( inum, jpdom_unknown, 'aisl'  , aisl )     
576         CALL iom_get( inum, jpdom_unknown, 'aislm1', aislm1 )     
577         CALL iom_close( inum )
578         ! Control print
579         IF(lwp) THEN
580            WRITE(numout,*)
581            WRITE(numout,*)' islmat: lecture aisl/aislm1 in numisl done'
582            WRITE(numout,*)' ~~~~~~'
583            WRITE(numout,*)
584            WRITE(numout,*) '        island matrix : '
585            WRITE(numout,*)
586           
587            DO jnj = 1, jpisl
588               WRITE(numout,'(8e12.4)') ( aisl(jni,jnj), jni = 1, jpisl )
589            END DO
590           
591            WRITE(numout,*)
592            WRITE(numout,*) '       inverse of the island matrix'
593            WRITE(numout,*)
594           
595            DO jnj = 1, jpisl
596               WRITE(numout,'(12e11.3)') ( aislm1(jni,jnj), jni=1,jpisl )
597            END DO
598         ENDIF
599 
600      ELSE
601 
602      ! II. Island matrix computation
603      ! =============================
604
605      DO jnj = 1, jpisl
606
607         ! Circulation of bsf(jnj) around island jni
608         
609         DO jj = 2, jpj
610               zwx(:,jj,1) = -( bsfisl(:,jj,jnj) - bsfisl(:,jj-1,jnj) )
611         END DO
612         zwx(:,1,1) = 0.e0
613         
614         DO jj = 1, jpj
615            DO ji = 2, jpi
616               zwx(ji,jj,2) = ( bsfisl(ji,jj,jnj) - bsfisl(ji-1,jj,jnj) )
617            END DO
618         END DO
619         zwx(1,:,2) = 0.e0
620         
621         ! Island matrix
622         
623         DO jni = 1, jpisl
624            aisl(jni,jnj) = 0.e0
625            DO jl = 1, 2
626               DO jj=1,jpj
627                  DO ji=1,jpi
628                     aisl(jni,jnj) = aisl(jni,jnj) + acisl2(ji,jj,jl,jni)*zwx(ji,jj,jl)
629                  END DO
630               END DO
631            END DO
632         END DO
633         
634      END DO
635      IF( lk_mpp ) THEN
636         DO jnj = 1, jpisl
637            DO jni = 1, jpisl
638               ztab(jni+(jnj-1)*jpisl) = aisl(jni,jnj)
639            END DO
640         END DO
641         CALL mpp_sum( ztab, jpisl*jpisl )   ! sum over the global domain
642      ENDIF
643
644      ! 1.3 Control print
645
646      IF(lwp) THEN
647         WRITE(numout,*)
648         WRITE(numout,*) 'islmat : island matrix'
649         WRITE(numout,*) '~~~~~~'
650         WRITE(numout,*)
651         
652         DO jnj = 1, jpisl
653            WRITE(numout,'(8e12.4)') ( aisl(jni,jnj), jni = 1, jpisl )
654         END DO
655      ENDIF
656     
657
658      ! 2. Invertion of the island matrix
659      ! ---------------------------------
660     
661      ! 2.1 Call of an imsl routine for the matrix invertion
662
663      CALL linrg( jpisl, aisl, jpisl, aislm1, jpisl )
664
665      ! 2.2 Control print
666     
667      IF(lwp) THEN
668         WRITE(numout,*)
669         WRITE(numout,*) 'islmat : inverse of the island matrix'
670         WRITE(numout,*) '~~~~~~'
671         WRITE(numout,*)
672         
673         DO jnj = 1, jpisl
674            WRITE(numout, '(12e11.3)')  '        ', ( aislm1(jni,jnj), jni=1, jpisl )
675         END DO
676      ENDIF
677     
678
679      ! 3. Output of aisl and aislm1 in numisl
680      ! --------------------------------------
681
682      CALL iom_rstput( 0, 0, inum, 'aisl'  , aisl   )
683      CALL iom_rstput( 0, 0, inum, 'aislm1', aislm1 )
684      CALL iom_close ( inum )
685
686      ENDIF
687
688   END SUBROUTINE isl_mat
689
690
691   SUBROUTINE isl_bsf
692      !!----------------------------------------------------------------------
693      !!                  ***  ROUTINE isl_bsf  ***
694      !!         
695      !! ** Purpose :
696      !!         Compute the barotropic stream function associated with each
697      !!      island using FETI or preconditioned conjugate gradient method
698      !!      or read them in numisl.
699      !!
700      !! ** Method :
701      !!
702      !! ** input/output file :
703      !!            numisl            : barotropic stream function associated
704      !!                          with each island of the domain
705      !!
706      !! ** Action :
707      !!       bsfisl, the streamfunction which takes the value 1 over island ni,
708      !!    and 0 over the others islands
709      !!       file 'numisl' barotropic stream function associated
710      !!                              with each island of the domain
711      !!
712      !! History :
713      !!        !  87-10  (G. Madec)  Original code
714      !!        !  91-11  (G. Madec)
715      !!        !  93-03  (G. Madec)  release 7.1
716      !!        !  93-04  (M. Guyon)  loops and suppress pointers
717      !!        !  96-11  (A. Weaver)  correction to preconditioning
718      !!        !  98-02  (M. Guyon)  FETI method
719      !!        !  99-11  (M. Imbard)  NetCDF FORMAT with IOIPSL
720      !!   8.5  !  02-08  (G. Madec)  Free form, F90
721      !!----------------------------------------------------------------------
722      !! * Modules used
723      USE ioipsl
724      USE iom
725      USE solpcg
726      USE solfet
727      USE solsor
728
729      !! * Local declarations
730      LOGICAL  ::   llog, llbon
731      CHARACTER (len=10) ::  clisl
732      CHARACTER (len=32) ::  clname = 'islands'
733      INTEGER  ::   &
734         inum                 ! temporary logical unit
735      INTEGER  ::   ji, jj, jni, jii, jnp  ! dummy loop indices
736      INTEGER  ::   iimlu, ijmlu, inmlu, iju
737      INTEGER  ::   ii, ij, icile, icut, inmax, indic
738      REAL(wp) ::   zepsr, zeplu, zgwgt
739      REAL(wp) ::   zep(jpisl), zprec(4)
740      REAL(wp) ::   zdt
741      REAL(wp) ::   t2p1(jpi,1,1)
742      INTEGER  ::   iloc
743      !!----------------------------------------------------------------------
744
745
746      ! 0. Initializations
747      ! ==================
748
749      icile         = 0      ! set to zero the convergence indicator
750      bsfisl(:,:,:) = 0.e0   ! set to zero of bsfisl
751
752
753      ! I. Lecture of bsfisl in numisl (if it exists)
754      ! =============================================
755
756      icut  = 0
757      iimlu = 0
758      ijmlu = 0
759      inmlu = 0
760      zeplu = 0.
761     
762      clname = 'islands'
763     
764      INQUIRE( FILE=clname, EXIST=llbon )
765! islands FILE does not EXIST : icut=999
766      IF( llbon ) THEN 
767
768         ! island FILE is present
769
770         CALL iom_open (clname, inum )
771         CALL iom_get( inum, jpdom_unknown, 'PRECISION', zprec )     
772
773         iimlu = NINT( zprec(1) )
774         ijmlu = NINT( zprec(2) )
775         inmlu = NINT( zprec(3) )
776         zeplu = zprec(4)
777         ! the read domain does not correspond to the model one : icut=999
778         IF( iimlu /= jpi .OR. ijmlu /= jpj .OR. inmlu /= jpisl ) THEN
779            icut = 999
780         ELSE
781            DO jni = 1, jpisl
782               IF( jni < 10 ) THEN
783                  WRITE(clisl,'("island",I1)') jni
784               ELSEIF( jni < 100 ) THEN
785                  WRITE(clisl,'("island",I2)') jni
786               ELSE
787                  WRITE(clisl,'("island",I3)') jni
788               ENDIF
789               CALL iom_get( inum, jpdom_autoglo, clisl, bsfisl(:,:,jni))     
790            END DO
791         ENDIF
792         CALL iom_close( inum )
793         !
794      ELSE 
795         ! islands FILE does not EXIST : icut=999
796         icut = 999
797      ENDIF
798   
799      ! the read precision is not the required one : icut=888
800      IF( zeplu > epsisl ) THEN
801         icut = 888
802      ENDIF
803
804      ! Control print
805      IF( icut == 999 ) THEN
806          IF(lwp) THEN
807              WRITE(numout,*)
808              WRITE(numout,*) 'islbsf : lecture bsfisl in numisl failed'
809              WRITE(numout,*) '~~~~~~'
810              WRITE(numout,*) '         icut= ', icut
811              WRITE(numout,*) '         imlu= ', iimlu, ' jmlu= ', ijmlu,   &
812                                      ' nilu= ', inmlu, ' epsisl lu= ', zeplu
813              WRITE(numout,*) '         the bsfisl are computed from zero'
814          ENDIF
815        ELSEIF( icut == 888 ) THEN
816          IF(lwp) THEN
817              WRITE(numout,*)
818              WRITE(numout,*) 'islbsf : lecture bsfisl in numisl done'
819              WRITE(numout,*) '~~~~~~'
820              WRITE(numout,*) '         the required accuracy is not reached'
821              WRITE(numout,*) '         epsisl lu= ', zeplu,' epsisl required ', epsisl
822              WRITE(numout,*) '         the bsfisl are computed from the read values'
823          ENDIF
824        ELSE
825          IF(lwp) THEN
826              WRITE(numout,*)
827              WRITE(numout,*) 'islbsf : lecture bsfisl in numisl done'
828              WRITE(numout,*) '~~~~~~'
829          ENDIF
830          RETURN
831      ENDIF
832     
833     
834      ! II. Compute the bsfisl (if icut=888 or 999)
835      ! ============================================
836
837      ! save nmax
838      inmax = nmax
839
840      ! set the number of iteration of island computation
841      nmax = nmisl
842
843      ! Loop over islands
844      ! -----------------
845     
846      DO jni = 1, jpisl
847
848
849         ! 1. Initalizations of island computation
850         ! ---------------------------------------
851         
852         ! Set the pcg solution gcb to zero
853         gcb(:,:) = 0.e0
854
855         ! Set first guess gcx either to zero or to the read bsfisl
856         IF( icut == 999 ) THEN
857            gcx(:,:) = 0.e0
858         ELSEIF( icut == 888 ) THEN
859            IF(lwp) WRITE(numout,*) ' islbsf: bsfisl read are used as first guess'
860            ! c a u t i o n: bsfisl masked because it contains 1 along island seaside
861            gcx(:,:) = bsfisl(:,:,jni) * bmask(:,:)
862         ENDIF
863         
864         ! Right hand side of the streamfunction equation
865         
866         IF( lk_mpp ) THEN
867
868            ! north fold treatment
869            IF( npolj == 3 .OR. npolj == 5)   iloc=jpiglo-(nimpp-1+nimppt(nono+1)-1)
870            IF( npolj == 4 .OR. npolj == 6)   iloc=jpiglo-2*(nimpp-1)
871            t2p1(:,1,1) = 0.e0
872            ! north and south grid-points
873            DO jii = 1, 2
874               DO jnp = 1, mnisl(jii,jni)
875                  ii = miisl(jnp,jii,jni)
876                  ij = mjisl(jnp,jii,jni)
877                  IF( ( npolj == 3 .OR. npolj == 4 ) .AND.   &
878                     ( ij == nlcj-1 .AND. jii == 1) ) THEN
879                     iju=iloc-ii+1
880                     t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
881                  ELSEIF( ( npolj == 5 .OR. npolj == 6 ) .AND.   &
882                     ( ij == nlcj-1 .AND. jii == 1) ) THEN
883                     iju=iloc-ii
884                     gcb(ii,ij) =  gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
885                     t2p1(iju,1,1) =  t2p1(iju,1,1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
886                  ELSE 
887                     gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
888                  ENDIF
889               END DO
890            END DO
891         
892            ! east and west grid-points
893         
894            DO jii = 3, 4
895               DO jnp = 1, mnisl(jii,jni)
896                  ii = miisl(jnp,jii,jni)
897                  ij = mjisl(jnp,jii,jni)
898                  gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
899               END DO
900            END DO
901
902            IF( lk_mpp )   CALL mpplnks( gcb )   !!bug ? should use an lbclnk ? is it possible?
903
904         ELSE
905
906            ! north and south grid-points
907            DO jii = 1, 2
908               DO jnp = 1, mnisl(jii,jni)
909                  ii = miisl(jnp,jii,jni)
910                  ij = mjisl(jnp,jii,jni)
911                  IF( ( nperio == 3 .OR. nperio == 4 ) .AND.   &
912                     ( ij == jpj-1 .AND. jii == 1) ) THEN
913                     gcb(jpi-ii+1,ij-1) = gcb(jpi-ii+1,ij-1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
914                  ELSEIF( ( nperio == 5 .OR. nperio == 6 ) .AND.   &
915                     ( ij == jpj-1 .AND. jii == 1) ) THEN
916                     gcb(ii,ij) =  gcb(ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
917                     gcb(jpi-ii,ij) = gcb(jpi-ii,ij) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij) 
918                  ELSE 
919                     gcb(ii,ij-jii+1) =  gcb(ii,ij-jii+1) + hur(ii,ij) * e1u(ii,ij) / e2u(ii,ij)
920                  ENDIF
921               END DO
922            END DO
923
924            ! east and west grid-points
925            DO jii = 3, 4
926               DO jnp = 1, mnisl(jii,jni)
927                  ii = miisl(jnp,jii,jni)
928                  ij = mjisl(jnp,jii,jni)
929                  IF( bmask(ii-jii+3,ij) /= 0. ) THEN
930                     gcb(ii-jii+3,ij) = gcb(ii-jii+3,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
931                  ELSE
932                     ! east-west cyclic boundary conditions
933                     IF( ii-jii+3 == 1 ) THEN
934                        gcb(jpim1,ij) = gcb(jpim1,ij) + hvr(ii,ij) * e2v(ii,ij) / e1v(ii,ij)
935                     ENDIF
936                  ENDIF
937               END DO
938            END DO
939         ENDIF
940
941         ! Preconditioned right hand side and absolute precision
942
943         IF( nsolv == 3 ) THEN 
944            ! FETI method
945            ncut    = 0
946            rnorme  = 0.e0
947            gcb(:,:) = bmask(:,:) * gcb(:,:)
948            DO jj = 1, jpj
949               DO ji = 1, jpi
950                  rnorme = rnorme + gcb(ji,jj) * gcb(ji,jj)
951               END DO
952            END DO
953           
954            IF( lk_mpp )   CALL mpp_sum( rnorme )
955
956            IF(lwp) WRITE(numout,*) 'rnorme ', rnorme
957            epsr  = epsisl * epsisl * rnorme
958            indic = 0
959         ELSE
960            ncut   = 0
961            rnorme = 0.e0
962            DO jj = 1, jpj
963               DO ji = 1, jpi
964                  gcb  (ji,jj) = gcdprc(ji,jj) * gcb(ji,jj)
965                  zgwgt  = gcdmat(ji,jj) * gcb(ji,jj)
966                  rnorme = rnorme + gcb(ji,jj) * zgwgt
967               END DO
968            END DO
969            IF( lk_mpp )   CALL mpp_sum( rnorme )   ! sum over the global domain
970
971            IF(lwp) WRITE(numout,*) 'rnorme ', rnorme
972            epsr = epsisl * epsisl * rnorme
973            indic = 0
974         ENDIF
975
976
977         ! 3. PCG solver for gcp.gcx=gcb in monotask
978         ! -------------------------------------------
979
980         IF( nsolv == 3 ) THEN
981            epsilo = epsisl       ! precision to compute Islands matrix A
982            CALL sol_fet( indic )  ! FETI method
983            epsilo = eps          ! precision to compute grad PS
984         ELSE
985            CALL sol_pcg( indic )  ! pcg method
986         ENDIF
987
988
989         ! 4. Save the solution in bsfisl
990         ! ------------------------------
991
992         bsfisl(:,:,jni) = gcx(:,:)
993
994
995         ! 5. Boundary conditions
996         ! ----------------------
997
998         ! set to 1. coastal gridpoints of the island
999         DO jnp = 1, mnisl(0,jni)
1000            ii = miisl(jnp,0,jni)
1001            ij = mjisl(jnp,0,jni)
1002            bsfisl(ii,ij,jni) = 1.
1003         END DO
1004
1005         ! cyclic boundary conditions
1006         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
1007            bsfisl( 1 ,:,jni) = bsfisl(jpim1,:,jni)
1008            bsfisl(jpi,:,jni) = bsfisl(  2  ,:,jni)
1009         ENDIF
1010         IF( nperio == 3 .OR. nperio == 4 ) THEN
1011            DO ji = 1, jpim1
1012               iju = jpi-ji+1
1013               bsfisl(ji,jpj-1,jni) = bsfisl(iju,jpj-2,jni)
1014               bsfisl(ji, jpj ,jni) = bsfisl(iju,jpj-3,jni)
1015            END DO
1016         ENDIF
1017         IF( nperio == 5 .OR. nperio == 6 ) THEN
1018            DO ji = 1, jpi-1
1019               iju=jpi-ji
1020               bsfisl(ji,jpj,jni) = bsfisl(iju,jpj-2,jni)
1021            END DO
1022            DO ji = jpi/2+1, jpi-1
1023               iju=jpi-ji
1024               bsfisl (ji,jpjm1,jni) = bsfisl(iju,jpjm1,jni)
1025            END DO
1026         ENDIF
1027         IF( lk_mpp )   CALL lbc_lnk( bsfisl(:,:,jni), 'G', 1. )   ! link at G-point
1028
1029
1030         ! 6. Control print
1031         ! ----------------
1032
1033         IF(lwp) WRITE(numout,*)
1034         IF(lwp) WRITE(numout,*) ' islbsf: island number: ', jni
1035         IF(lwp) WRITE (numout,9290) niter, res, SQRT(epsr)/epsisl
1036         zep(jni) =  MAX(epsisl, res/(SQRT(epsr)/epsisl))
1037         IF( indic <  0 ) THEN
1038            icile = icile-1
1039            IF(lwp) WRITE(numout,*) '  pcg do not converge for island: ', jni
1040            IF(lwp) WRITE(numout,*) '      Precision reached: ',zep(jni)
1041         ENDIF
1042
10439290     FORMAT('        niter :',i4,' , res :',e20.10,' , gcb :',e20.10)
1044
1045         !                                          !====================
1046      END DO                                        !  End Loop islands
1047      !                                             !====================
1048
1049      ! 7. Reset PCG
1050      ! ------------
1051
1052      ! reset the number of iteration for pcg
1053      nmax = inmax
1054
1055      ! reset to zero pcg arrays
1056      gcx  (:,:) = 0.e0
1057      gcxb (:,:) = 0.e0
1058      gcb  (:,:) = 0.e0
1059      gcr  (:,:) = 0.e0
1060      gcdes(:,:) = 0.e0
1061      gccd (:,:) = 0.e0
1062
1063
1064      ! III. Output of bsfisl in numisl
1065      ! ===============================
1066
1067      zprec(1) = FLOAT(jpi)
1068      zprec(2) = FLOAT(jpj)
1069      zprec(3) = FLOAT(jpisl)
1070      IF(lwp) WRITE(numout,*) clname
1071      CALL iom_open  ( clname, numisl, ldwrt = .TRUE., kiolib = jprstlib )
1072      IF( icile == 0 .AND. icut /= 0 ) THEN
1073         IF(lwp) THEN
1074            WRITE(numout,*)
1075            WRITE(numout,*)' islbsf: write bsfisl in numisl ', numisl
1076            WRITE(numout,*)' -------------------------------------'
1077         ENDIF
1078         zprec(4) = epsisl
1079         CALL iom_rstput( 0, 0, numisl, 'PRECISION' , zprec )
1080         DO jni = 1, jpisl
1081            IF(jni < 10) THEN
1082               WRITE(clisl,'("island",I1)') jni
1083            ELSE IF(jni < 100) THEN
1084               WRITE(clisl,'("island",I2)') jni
1085            ELSE
1086               WRITE(clisl,'("island",I3)') jni
1087            ENDIF
1088            CALL iom_rstput( 0, 0, numisl, clisl  , bsfisl(:,:,jni) )
1089         END DO
1090         CALL iom_close( numisl )
1091      ENDIF
1092
1093      IF( icile < 0 ) THEN
1094         IF(lwp) THEN
1095            WRITE(numout,*)
1096            WRITE(numout,*) ' islbsf: number of island without convergence : ',ABS(icile)
1097            WRITE(numout,*) ' ---------------------------------------------'
1098         ENDIF
1099         zepsr = epsisl
1100         DO jni = 1, jpisl
1101            IF(lwp) WRITE(numout,*) '    isl ',jni,' precision reached ', zep(jni)
1102            zepsr = MAX( zep(jni), zepsr )
1103         END DO
1104         IF( zepsr == 0. ) zepsr = epsisl
1105         IF(lwp) THEN
1106            WRITE(numout,*) ' save value of precision reached: ',zepsr
1107            WRITE(numout,*)
1108            WRITE(numout,*)' islbsf: save bsfisl in numisl ',numisl
1109            WRITE(numout,*)' -------------------------------------'
1110         ENDIF
1111
1112         zprec(4) = zepsr
1113         CALL iom_rstput( 0, 0, numisl, 'PRECISION' , zprec )
1114         DO jni = 1, jpisl
1115            IF( jni < 10 ) THEN
1116               WRITE(clisl,'("island",I1)') jni
1117            ELSE IF( jni < 100 ) THEN
1118               WRITE(clisl,'("island",I2)') jni
1119            ELSE
1120               WRITE(clisl,'("island",I3)') jni
1121            ENDIF
1122            CALL iom_rstput( 0, 0, numisl, clisl  , bsfisl(:,:,jni) )
1123         END DO
1124         CALL iom_close( numisl )
1125         CALL ctl_stop( ' ' )
1126      ENDIF
1127
1128   END SUBROUTINE isl_bsf
1129
1130
1131   SUBROUTINE isl_dyn_spg
1132      !!----------------------------------------------------------------------
1133      !!                  ***  routine isl_dyn_spg  ***
1134      !!
1135      !! ** Purpose :   Compute and add the island contribution to the
1136      !!      barotropic stream function trend.
1137      !!
1138      !!
1139      !! ** Method  :   Rigid-lid appromimation: ...????
1140      !!
1141      !! ** Action : - Update bsfd with the island contribution
1142      !!
1143      !! History :
1144      !!   9.0  !  03-09  (G. Madec)  isolate island computation
1145      !!---------------------------------------------------------------------
1146
1147      !! * Local declarations
1148      INTEGER ::   ji, jj, jni, jnj    ! dummy loop indices
1149      !!----------------------------------------------------------------------
1150     
1151
1152      ! compute the island potential
1153      ! ----------------------------
1154      DO jni = 1, jpisl                      ! second member
1155         bisl(jni) =  0.e0
1156         DO jj = 2, jpj
1157            DO ji = 2, jpi
1158               bisl(jni) =  bisl(jni) + acisl1(ji,jj,1,jni) *  spgu(ji,jj)                  &
1159                  &                   + acisl1(ji,jj,2,jni) *  spgv(ji,jj)                  &
1160                  &                   + acisl2(ji,jj,1,jni) * ( gcx(ji,jj)-gcx(ji,jj-1) )   &
1161                  &                   - acisl2(ji,jj,2,jni) * ( gcx(ji,jj)-gcx(ji-1,jj) )
1162            END DO
1163         END DO
1164      END DO
1165      IF( lk_mpp )   CALL mpp_sum( bisl, jpisl )   ! sum over the global domain
1166
1167      DO jni = 1, jpisl                     ! Island stream function trend
1168         visl(jni) = 0.e0
1169         DO jnj = 1, jpisl
1170            visl(jni) = visl(jni) + aislm1(jni,jnj) * bisl(jnj)
1171         END DO
1172      END DO
1173
1174      ! update the bsf trend ( caution : bsfd is not zero along island coastlines, dont mask it ! )
1175      ! --------------------
1176      DO jj = 1, jpj
1177         DO jni = 1, jpisl
1178            DO ji = 1, jpi
1179               bsfd(ji,jj) = bsfd(ji,jj) + visl(jni) * bsfisl(ji,jj,jni)
1180            END DO
1181         END DO
1182      END DO
1183
1184   END SUBROUTINE isl_dyn_spg
1185
1186
1187   SUBROUTINE isl_stp_ctl( kt, kindic )
1188      !!----------------------------------------------------------------------
1189      !!                    ***  ROUTINE isl_stp_ctl  ***
1190      !!                     
1191      !! ** Purpose :   ???
1192      !!
1193      !! ** Method  : - print island potential
1194      !!
1195      !! History :
1196      !!   9.0  !  03-09  (G. Madec)  isolated from stp_ctl
1197      !!----------------------------------------------------------------------
1198      !! * Arguments
1199      INTEGER, INTENT(in   ) ::   kt        ! ocean time-step index
1200      INTEGER, INTENT(inout) ::   kindic    ! indicator of solver convergence
1201
1202      !! * local declarations
1203      INTEGER  ::   jni                     ! dummy loop indice
1204      REAL(wp) ::   zfact                   ! temporary scalar
1205      !!----------------------------------------------------------------------
1206      !!  OPA 8.5, LODYC-IPSL (2002)
1207      !!----------------------------------------------------------------------
1208
1209      IF( kt == nit000 .AND. lwp ) THEN
1210         WRITE(numout,*)
1211         WRITE(numout,*) 'isl_stp_ctl : time-stepping control'
1212         WRITE(numout,*) '~~~~~~~~~~~'
1213      ENDIF
1214
1215      ! Island trends
1216      DO jni = 1, jpisl
1217         zfact = 0.
1218         IF( miisl(1,0,jni) /= 0 .AND. mjisl(1,0,jni) /= 0 ) THEN
1219            zfact = 1.e-6 * bsfn(miisl(1,0,jni),mjisl(1,0,jni))
1220         ENDIF
1221         IF( lk_mpp )   CALL mpp_isl( zfact )
1222
1223         IF(lwp) WRITE(numisp,9300) kt, jni, zfact, visl(jni)
1224         IF( MOD( kt, nwrite ) == 0 .OR. kindic < 0     &
1225            .OR. ( kt == nit000 .AND. kindic > 0 )      &
1226            .OR.   kt == nitend                         ) THEN
1227            IF( jni == 1 .AND. lwp ) THEN
1228               WRITE(numout,*)
1229               WRITE(numout,*) 'isl_stp_ctl : island bsf'
1230               WRITE(numout,*) '~~~~~~~~~~~'
1231            ENDIF
1232            IF(lwp) WRITE(numout,9300) kt, jni, zfact, visl(jni)
1233         ENDIF
1234      END DO
12359300  FORMAT(' it : ',i8,' island  :',i4,'   BSF (Sverdrup) : ',f7.2, '   visl : ',e15.6)
1236
1237   END SUBROUTINE isl_stp_ctl
1238
1239#else
1240   !!----------------------------------------------------------------------
1241   !!   Default option                                         Empty module
1242   !!----------------------------------------------------------------------
1243   LOGICAL, PUBLIC, PARAMETER ::   lk_isl = .FALSE.    !: 'key_islands' flag
1244CONTAINS
1245   SUBROUTINE isl_dom                        ! Empty routine
1246   END SUBROUTINE isl_dom
1247   SUBROUTINE isl_bsf                        ! Empty routine
1248   END SUBROUTINE isl_bsf
1249   SUBROUTINE isl_mat                        ! Empty routine
1250   END SUBROUTINE isl_mat
1251   SUBROUTINE isl_dyn_spg                    ! Empty routine
1252   END SUBROUTINE isl_dyn_spg
1253   SUBROUTINE isl_stp_ctl( kt, kindic )      ! Empty routine
1254      WRITE(*,*) 'isl_stp_ctl: You should not have seen this print! error?', kt, kindic
1255   END SUBROUTINE isl_stp_ctl
1256#endif
1257
1258   !!======================================================================
1259END MODULE solisl
Note: See TracBrowser for help on using the repository browser.