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

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

Initial revision

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