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

source: tags/nemo_dev_x6/NEMO/OPA_SRC/SOL/solisl.F90 @ 3532

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

CT : BUGFIX028 : # Add USE modules and correct calling sequences CALL mmp_sum() and CALL sol_pcg() CALL sol_fet()

# Use logical key "lk_isl" instead of "l_isl"

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