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

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

CT : UPDATE001 : First major NEMO update

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