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

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

CT : UPDATE172 : remove all direct acces modules and the related cpp key key_fdir

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