source: trunk/NEMO/LIM_SRC/limrhg.F90 @ 12

Last change on this file since 12 was 12, checked in by opalod, 17 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: 35.8 KB
Line 
1MODULE limrhg
2#if defined key_ice_lim
3   !!======================================================================
4   !!                     ***  MODULE  limrhg  ***
5   !!   Ice rheology :  performs sea ice rheology
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   lim_rhg   : computes ice velocities
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE phycst
13   USE ice_oce         ! ice variables
14   USE dom_ice
15   USE ice
16   USE lbclnk
17   USE in_out_manager  ! I/O manager
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Routine accessibility
23   PUBLIC lim_rhg  ! routine called by lim_dyn
24
25   !! * Module variables
26     REAL(wp)  ::            &  ! constant values
27         rzero   = 0.e0   ,  &
28         rone    = 1.e0
29   !!----------------------------------------------------------------------
30   !!   LIM 2.0 , UCL-LODYC-IPSL  (2003)
31   !!----------------------------------------------------------------------
32
33CONTAINS
34
35   SUBROUTINE lim_rhg( khemi )
36      !!-------------------------------------------------------------------
37      !!                 ***  SUBROUTINR lim_rhg  ***
38      !! ** purpose :   determines the velocity field of sea ice by using
39      !!  atmospheric (wind stress) and oceanic (water stress and surface
40      !!  tilt) forcings. Ice-ice interaction is described by a non-linear
41      !!  viscous-plastic law including shear strength and a bulk rheology.   
42      !!
43      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity
44      !!
45      !! History :
46      !!   0.0  !  93-12  (M.A. Morales Maqueda.)  Original code
47      !!   1.0  !  94-12  (H. Goosse)
48      !!   2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp
49      !!-------------------------------------------------------------------
50      ! * Arguments
51      INTEGER, INTENT(in) ::   khemi   ! -1/1 = South/North hemisphere flag
52
53      ! * Local variables
54      INTEGER ::   ji, jj              ! dummy loop indices
55
56      INTEGER  :: &
57         i_j1, i_j2, i_jpj, i_jpjm1, & ! ????
58         iim1, ijm1, iip1 , ijp1   , & ! temporary integers
59         iter, jter                    !    "          "
60
61      REAL(wp) :: &
62         ze11  , ze12  , ze22  , ze21  ,   &  ! temporary scalars
63         zt11  , zt12  , zt21  , zt22  ,   &  !    "         "
64         zvis11, zvis21, zvis12, zvis22,   &  !    "         "
65         zgphsx, ztagnx, zusw  ,           &  !    "         "
66         zgphsy, ztagny                       !    "         "
67      REAL(wp) :: &
68         zresm, zunw, zvnw, zur, zvr, zmod, za, zac, &
69         zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal,  &
70         ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag
71      REAL(wp),DIMENSION(jpj) ::   &
72         zind                                 ! i-averaged indicator of sea-ice
73      REAL(wp),DIMENSION(jpi,jpj) ::   &
74         zpresh, zfrld, zmass, zcorl,     &
75         zu0, zv0, zviszeta, zviseta,     &
76         zc1u, zc1v, zc2u, zc2v, za1ct, za2ct, za1, za2, zb1, zb2,  &
77         zc1, zc2, zd1, zd2, zden, zu_ice, zv_ice, zresr
78      REAL(wp),DIMENSION(jpi,jpj,2,2) :: &
79         zsigm11, zsigm12, zsigm22, zsigm21
80      !!-------------------------------------------------------------------
81     
82      !  Store initial velocities.
83      zu0(:,:) = u_ice(:,:)
84      zv0(:,:) = v_ice(:,:)
85
86      ! Numerical characteristics.                                       
87      ! --------------------------
88
89      !  Define the j-limits where ice dynamics is computed
90      ! ---------------------------------------------------
91
92      DO jj = 1, jpj
93         zind(jj) = SUM( frld(:,jj) )   ! = FLOAT(jpj) if ocean everwhere on a j-line
94      END DO
95
96      IF( khemi ==  1 )  THEN      ! Northern hemisphere
97         i_j1  = jeq
98         i_jpj = jpj
99         DO jj = jpj, jeq, -1
100            IF( zind(jj) < FLOAT(jpi) )   i_j1 = jj
101         END DO
102         i_j1 = MAX( 1, i_j1-1)
103         IF( l_ctl .AND. lwp ) WRITE(numout,*) 'lim_rhg : NH i_j1 = ', i_j1, ' ij_pj = ', i_jpj
104      ELSE                         ! Southern hemisphere
105         i_j1  = 2
106         i_jpj = jpj
107         DO jj = 1, jeq
108            IF( zind(jj) < FLOAT(jpi) ) i_jpj = jj
109         END DO
110         i_jpj = MIN( jpj, i_jpj+2)
111         IF( l_ctl .AND. lwp ) WRITE(numout,*) 'lim_rhg : SH i_j1 = ', i_j1, ' ij_pj = ', i_jpj
112      ENDIF
113      i_j2    = i_j1  + 1   ! inner domain indices
114      i_jpjm1 = i_jpj - 1
115
116
117      !  2) Sign of turning angle for oceanic drag.                          |
118      !-----------------------------------------------------------------------
119
120      zsang = REAL( khemi ) * sangvg   ! only the sinus changes its sign with the hemisphere
121
122
123      !  3) Ice mass, ice strength, and wind stress at the center            |
124      !     of the grid squares.                                             |
125      !-----------------------------------------------------------------------
126
127      DO jj = i_j1 , i_jpjm1
128         DO ji = 1 , jpi
129            za1(ji,jj)    = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
130            zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
131#if defined key_lim_cp1 && defined key_coupled
132            zb1(ji,jj)    = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) )
133            zb2(ji,jj)    = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) )
134#else
135            zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
136            zb2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
137#endif
138         END DO
139      END DO
140
141
142      !---------------------------------------------------------------------------
143      !  Wind stress, coriolis and mass terms at the corners of the grid squares |
144      !  Gradient of ice strenght.                                               |
145      !---------------------------------------------------------------------------
146         
147      DO jj = i_j2, i_jpjm1
148         DO ji = 2, jpi
149            zstms = tms(ji,jj  ) * wght(ji,jj,2,2) + tms(ji-1,jj  ) * wght(ji,jj,1,2)   &
150               &  + tms(ji,jj-1) * wght(ji,jj,2,1) + tms(ji-1,jj-1) * wght(ji,jj,1,1)
151            zusw  = 1.0 / MAX( zstms, epsd )
152
153            zt11 = tms(ji  ,jj  ) * frld(ji  ,jj  ) 
154            zt12 = tms(ji-1,jj  ) * frld(ji-1,jj  ) 
155            zt21 = tms(ji  ,jj-1) * frld(ji  ,jj-1) 
156            zt22 = tms(ji-1,jj-1) * frld(ji-1,jj-1)
157
158            ! Leads area.
159            zfrld(ji,jj) =  (  zt11 * wght(ji,jj,2,2) + zt12 * wght(ji,jj,1,2)   &
160               &             + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw
161
162            ! Mass and coriolis coeff.
163            zmass(ji,jj) = ( za1(ji,jj  ) * wght(ji,jj,2,2) + za1(ji-1,jj  ) * wght(ji,jj,1,2)   &
164               &           + za1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
165            zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj)
166
167            ! Wind stress.
168#if defined key_lim_cp1 && defined key_coupled
169            ztagnx = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   &
170               &     + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
171            ztagny = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   &
172               &     + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
173#else
174            ztagnx = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   &
175               &     + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj)
176            ztagny = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   &
177               &     + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj)
178#endif
179
180            ! Gradient of ice strength
181            zgphsx =   ( alambd(ji,jj,2,2,2,1) - alambd(ji,jj,2,1,2,1) ) * zpresh(ji  ,jj-1)   &
182               &     + ( alambd(ji,jj,2,2,2,2) - alambd(ji,jj,2,1,2,2) ) * zpresh(ji  ,jj  )   &
183               &     - ( alambd(ji,jj,2,2,1,1) + alambd(ji,jj,2,1,1,1) ) * zpresh(ji-1,jj-1)   &
184               &     - ( alambd(ji,jj,2,2,1,2) + alambd(ji,jj,2,1,1,2) ) * zpresh(ji-1,jj  )
185
186            zgphsy = - ( alambd(ji,jj,1,1,2,1) + alambd(ji,jj,1,2,2,1) ) * zpresh(ji  ,jj-1)   &
187               &     - ( alambd(ji,jj,1,1,1,1) + alambd(ji,jj,1,2,1,1) ) * zpresh(ji-1,jj-1)   &
188               &     + ( alambd(ji,jj,1,1,2,2) - alambd(ji,jj,1,2,2,2) ) * zpresh(ji  ,jj  )   &
189               &     + ( alambd(ji,jj,1,1,1,2) - alambd(ji,jj,1,2,1,2) ) * zpresh(ji-1,jj  )
190
191            ! Computation of the velocity field taking into account the ice-ice interaction.                                 
192            ! Terms that are independent of the velocity field.
193            za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_oce(ji,jj) - zgphsx
194            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_oce(ji,jj) - zgphsy
195         END DO
196      END DO
197
198!! inutile!!
199!!??    CALL lbc_lnk( za1ct, 'I', -1. )
200!!??    CALL lbc_lnk( za2ct, 'I', -1. )
201
202
203      ! SOLUTION OF THE MOMENTUM EQUATION.
204      !------------------------------------------
205      !                                                   ! ==================== !
206      DO iter = 1 , 2 * nbiter                            !    loop over iter    !
207         !                                                ! ==================== !       
208         zindu = MOD( iter , 2 )
209         zusdtp = ( zindu * 2.0 + ( 1.0 - zindu ) * 1.0 )  * REAL( nbiter ) / rdt_ice
210
211         ! Computation of free drift field for free slip boundary conditions.
212
213           DO jj = i_j1, i_jpjm1
214              DO ji = 1, jpim1
215                 !- Rate of strain tensor.
216                 zt11 =   akappa(ji,jj,1,1) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) - u_ice(ji,jj  ) - u_ice(ji  ,jj+1) )  &
217                    &   + akappa(ji,jj,1,2) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) + v_ice(ji,jj  ) + v_ice(ji  ,jj+1) )
218                 zt12 = - akappa(ji,jj,2,2) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) - u_ice(ji,jj+1) - u_ice(ji+1,jj+1) )  &
219                    &   - akappa(ji,jj,2,1) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) + v_ice(ji,jj+1) + v_ice(ji+1,jj+1) )
220                 zt22 = - akappa(ji,jj,2,2) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) - v_ice(ji,jj+1) - v_ice(ji+1,jj+1) )  &
221                    &   + akappa(ji,jj,2,1) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) + u_ice(ji,jj+1) + u_ice(ji+1,jj+1) )
222                 zt21 =   akappa(ji,jj,1,1) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) - v_ice(ji,jj  ) - v_ice(ji  ,jj+1) )  &
223                    &   - akappa(ji,jj,1,2) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) + u_ice(ji,jj  ) + u_ice(ji  ,jj+1) )
224
225                 !- Rate of strain tensor.
226                 zdgp = zt11 + zt22
227                 zdgi = zt12 + zt21
228                 ztrace2 = zdgp * zdgp 
229                 zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi
230
231                 !  Creep limit depends on the size of the grid.
232                 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2),  creepl)
233
234                 !-  Computation of viscosities.
235                 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn )
236                 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2
237              END DO
238           END DO
239!!??       CALL lbc_lnk( zviszeta, 'I', -1. )  ! or T point???   semble reellement inutile
240!!??       CALL lbc_lnk( zviseta , 'I', -1. )
241
242
243           !-  Determination of zc1u, zc2u, zc1v and zc2v.
244           DO jj = i_j2, i_jpjm1
245              DO ji = 2, jpim1
246                 ze11   =  akappa(ji-1,jj-1,1,1)
247                 ze12   = +akappa(ji-1,jj-1,2,2)
248                 ze22   =  akappa(ji-1,jj-1,2,1)
249                 ze21   = -akappa(ji-1,jj-1,1,2)
250                 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
251                 zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
252                 zvis12 =       zviseta (ji-1,jj-1) + dm
253                 zvis21 =       zviseta (ji-1,jj-1)
254
255                 zdiag = zvis22 * ( ze11 + ze22 )
256                 zsigm11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
257                 zsigm12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
258                 zsigm22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
259                 zsigm21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
260
261                 ze11   = -akappa(ji,jj-1,1,1)
262                 ze12   = +akappa(ji,jj-1,2,2)
263                 ze22   =  akappa(ji,jj-1,2,1)
264                 ze21   = -akappa(ji,jj-1,1,2)
265                 zvis11 = 2.0 * zviseta (ji,jj-1) + dm
266                 zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
267                 zvis12 =       zviseta (ji,jj-1) + dm
268                 zvis21 =       zviseta (ji,jj-1)
269
270                 zdiag = zvis22 * ( ze11 + ze22 )
271                 zsigm11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
272                 zsigm12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
273                 zsigm22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
274                 zsigm21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
275
276                 ze11   =  akappa(ji-1,jj,1,1)
277                 ze12   = -akappa(ji-1,jj,2,2)
278                 ze22   =  akappa(ji-1,jj,2,1)
279                 ze21   = -akappa(ji-1,jj,1,2)
280                 zvis11 = 2.0 * zviseta (ji-1,jj) + dm
281                 zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
282                 zvis12 =       zviseta (ji-1,jj) + dm
283                 zvis21 =       zviseta (ji-1,jj)
284
285                 zdiag = zvis22 * ( ze11 + ze22 ) 
286                 zsigm11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
287                 zsigm12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
288                 zsigm22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
289                 zsigm21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
290
291                 ze11   = -akappa(ji,jj,1,1)
292                 ze12   = -akappa(ji,jj,2,2)
293                 ze22   =  akappa(ji,jj,2,1)
294                 ze21   = -akappa(ji,jj,1,2)
295                 zvis11 = 2.0 * zviseta (ji,jj) + dm
296                 zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
297                 zvis12 =       zviseta (ji,jj) + dm
298                 zvis21 =       zviseta (ji,jj)
299
300                 zdiag = zvis22 * ( ze11 + ze22 )
301                 zsigm11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
302                 zsigm12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
303                 zsigm22(ji,jj,2,2) =  zvis11 * ze22 + zdiag 
304                 zsigm21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
305              END DO
306           END DO
307
308           DO jj = i_j2, i_jpjm1
309              DO ji = 2, jpim1
310                 zc1u(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm11(ji,jj,2,2)   &
311                    &        - alambd(ji,jj,2,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm11(ji,jj,1,2)   &
312                    &        - alambd(ji,jj,1,1,2,1) * zsigm12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm12(ji,jj,1,1)   &
313                    &        + alambd(ji,jj,1,1,2,2) * zsigm12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm12(ji,jj,1,2)   &
314                    &        + alambd(ji,jj,1,2,1,1) * zsigm21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zsigm21(ji,jj,2,1)   &
315                    &        + alambd(ji,jj,1,2,1,2) * zsigm21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zsigm21(ji,jj,2,2)   &
316                    &        - alambd(ji,jj,2,1,1,1) * zsigm22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zsigm22(ji,jj,2,1)   &
317                    &        - alambd(ji,jj,2,1,1,2) * zsigm22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zsigm22(ji,jj,2,2)
318                 
319                 zc2u(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm21(ji,jj,2,2)   &
320                    &        - alambd(ji,jj,2,2,1,1) * zsigm21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm21(ji,jj,1,2)   &
321                    &        - alambd(ji,jj,1,1,2,1) * zsigm22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm22(ji,jj,1,1)   &
322                    &        + alambd(ji,jj,1,1,2,2) * zsigm22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm22(ji,jj,1,2)   &
323                    &        - alambd(ji,jj,1,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zsigm11(ji,jj,2,1)   &
324                    &        - alambd(ji,jj,1,2,1,2) * zsigm11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zsigm11(ji,jj,2,2)   &
325                    &        + alambd(ji,jj,2,1,1,1) * zsigm12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zsigm12(ji,jj,2,1)   &
326                    &        + alambd(ji,jj,2,1,1,2) * zsigm12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zsigm12(ji,jj,2,2)
327             END DO
328           END DO
329
330           DO jj = i_j2, i_jpjm1
331              DO ji = 2, jpim1
332                 !  zc1v , zc2v.
333                 ze11   =  akappa(ji-1,jj-1,1,2)
334                 ze12   = -akappa(ji-1,jj-1,2,1)
335                 ze22   = +akappa(ji-1,jj-1,2,2)
336                 ze21   =  akappa(ji-1,jj-1,1,1)
337                 zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
338                 zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
339                 zvis12 =       zviseta (ji-1,jj-1) + dm
340                 zvis21 =       zviseta (ji-1,jj-1)
341
342                 zdiag = zvis22 * ( ze11 + ze22 )
343                 zsigm11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
344                 zsigm12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
345                 zsigm22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
346                 zsigm21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
347 
348                 ze11   =  akappa(ji,jj-1,1,2)
349                 ze12   = -akappa(ji,jj-1,2,1)
350                 ze22   = +akappa(ji,jj-1,2,2)
351                 ze21   = -akappa(ji,jj-1,1,1)
352                 zvis11 = 2.0 * zviseta (ji,jj-1) + dm
353                 zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
354                 zvis12 =       zviseta (ji,jj-1) + dm
355                 zvis21 =       zviseta (ji,jj-1)
356
357                 zdiag = zvis22 * ( ze11 + ze22 )
358                 zsigm11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
359                 zsigm12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
360                 zsigm22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
361                 zsigm21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
362
363                 ze11   =  akappa(ji-1,jj,1,2)
364                 ze12   = -akappa(ji-1,jj,2,1)
365                 ze22   = -akappa(ji-1,jj,2,2)
366                 ze21   =  akappa(ji-1,jj,1,1)
367                 zvis11 = 2.0 * zviseta (ji-1,jj) + dm
368                 zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
369                 zvis12 =       zviseta (ji-1,jj) + dm
370                 zvis21 =       zviseta (ji-1,jj)
371
372                 zdiag = zvis22 * ( ze11 + ze22 )
373                 zsigm11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
374                 zsigm12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
375                 zsigm22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
376                 zsigm21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
377
378                 ze11   =  akappa(ji,jj,1,2)
379                 ze12   = -akappa(ji,jj,2,1)
380                 ze22   = -akappa(ji,jj,2,2)
381                 ze21   = -akappa(ji,jj,1,1)
382                 zvis11 = 2.0 * zviseta (ji,jj) + dm
383                 zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
384                 zvis12 =       zviseta (ji,jj) + dm
385                 zvis21 =       zviseta (ji,jj)
386
387                 zdiag = zvis22 * ( ze11 + ze22 )
388                 zsigm11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
389                 zsigm12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
390                 zsigm22(ji,jj,2,2) =  zvis11 * ze22 + zdiag
391                 zsigm21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
392
393              END DO
394           END DO
395
396           DO jj = i_j2, i_jpjm1
397              DO ji = 2, jpim1
398                 zc1v(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm11(ji,jj,2,2)   &
399                    &        - alambd(ji,jj,2,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm11(ji,jj,1,2)   &
400                    &        - alambd(ji,jj,1,1,2,1) * zsigm12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm12(ji,jj,1,1)   &
401                    &        + alambd(ji,jj,1,1,2,2) * zsigm12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm12(ji,jj,1,2)   &
402                    &        + alambd(ji,jj,1,2,1,1) * zsigm21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zsigm21(ji,jj,2,1)   &
403                    &        + alambd(ji,jj,1,2,1,2) * zsigm21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zsigm21(ji,jj,2,2)   &
404                    &        - alambd(ji,jj,2,1,1,1) * zsigm22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zsigm22(ji,jj,2,1)   &
405                    &        - alambd(ji,jj,2,1,1,2) * zsigm22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zsigm22(ji,jj,2,2)
406                 zc2v(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm21(ji,jj,2,2)   &
407                    &        - alambd(ji,jj,2,2,1,1) * zsigm21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm21(ji,jj,1,2)   &
408                    &        - alambd(ji,jj,1,1,2,1) * zsigm22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm22(ji,jj,1,1)   &
409                    &        + alambd(ji,jj,1,1,2,2) * zsigm22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm22(ji,jj,1,2)   &
410                    &        - alambd(ji,jj,1,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zsigm11(ji,jj,2,1)   &
411                    &        - alambd(ji,jj,1,2,1,2) * zsigm11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zsigm11(ji,jj,2,2)   &
412                    &        + alambd(ji,jj,2,1,1,1) * zsigm12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zsigm12(ji,jj,2,1)   &
413                    &        + alambd(ji,jj,2,1,1,2) * zsigm12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zsigm12(ji,jj,2,2)
414              END DO
415           END DO
416
417         ! Relaxation.
418           
419iflag:   DO jter = 1 , nbitdr
420
421            !  Store previous drift field.   
422            DO jj = i_j1, i_jpjm1
423               zu_ice(:,jj) = u_ice(:,jj)
424               zv_ice(:,jj) = v_ice(:,jj)
425            END DO
426
427            DO jj = i_j2, i_jpjm1
428                DO ji = 2, jpim1
429                  zur     = u_ice(ji,jj) - u_oce(ji,jj)
430                  zvr     = v_ice(ji,jj) - v_oce(ji,jj)
431                  zmod    = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) )
432                  za      = rhoco * zmod
433                  zac     = za * cangvg
434                  zmpzas  = alpha * zcorl(ji,jj) + za * zsang
435                  zmassdt = zusdtp * zmass(ji,jj)
436                  zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj)
437
438                  za1(ji,jj) =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   &
439                     &        + za * ( cangvg * u_oce(ji,jj) - zsang * v_oce(ji,jj) )
440
441                  za2(ji,jj) =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   &
442                     &        + za * ( cangvg * v_oce(ji,jj) + zsang * u_oce(ji,jj) )
443
444                  zb1(ji,jj)  = zmassdt + zac - zc1u(ji,jj)
445                  zb2(ji,jj)  = zmpzas  - zc2u(ji,jj)
446                  zc1(ji,jj)  = zmpzas  + zc1v(ji,jj)
447                  zc2(ji,jj)  = zmassdt + zac  - zc2v(ji,jj) 
448                  zdeter      = zc1(ji,jj) * zb2(ji,jj) + zc2(ji,jj) * zb1(ji,jj)
449                  zden(ji,jj) = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) )
450               END DO
451            END DO
452
453            ! The computation of ice interaction term is splitted into two parts
454            !-------------------------------------------------------------------------
455
456            ! Terms that do not involve already up-dated velocities.
457         
458            DO jj = i_j2, i_jpjm1
459               DO ji = 2, jpim1
460                  iim1 = ji
461                  ijm1 = jj - 1
462                  iip1 = ji + 1
463                  ijp1 = jj
464                  ze11 =   akappa(iim1,ijm1,1,1) * u_ice(iip1,ijp1) + akappa(iim1,ijm1,1,2) * v_ice(iip1,ijp1)
465                  ze12 = + akappa(iim1,ijm1,2,2) * u_ice(iip1,ijp1) - akappa(iim1,ijm1,2,1) * v_ice(iip1,ijp1)
466                  ze22 = + akappa(iim1,ijm1,2,2) * v_ice(iip1,ijp1) + akappa(iim1,ijm1,2,1) * u_ice(iip1,ijp1)
467                  ze21 =   akappa(iim1,ijm1,1,1) * v_ice(iip1,ijp1) - akappa(iim1,ijm1,1,2) * u_ice(iip1,ijp1)
468                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
469                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
470                  zvis12 =       zviseta (iim1,ijm1) + dm
471                  zvis21 =       zviseta (iim1,ijm1)
472                  zdiag = zvis22 * ( ze11 + ze22 )
473                  zsigm11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
474                  zsigm12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
475                  zsigm22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
476                  zsigm21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
477
478
479                  iim1 = ji - 1
480                  ijm1 = jj
481                  iip1 = ji
482                  ijp1 = jj + 1                   
483                  ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   &
484                     &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijp1) + v_ice(iim1,ijp1) )
485                  ze12 = + akappa(iim1,ijm1,2,2) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )   &
486                     &   - akappa(iim1,ijm1,2,1) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )
487                  ze22 = + akappa(iim1,ijm1,2,2) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )   &
488                     &   + akappa(iim1,ijm1,2,1) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )
489                  ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   &
490                     &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijp1) + u_ice(iim1,ijp1) )
491                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
492                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
493                  zvis12 =       zviseta (iim1,ijm1) + dm
494                  zvis21 =       zviseta (iim1,ijm1)
495                  zdiag = zvis22 * ( ze11 + ze22 )
496                  zsigm11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
497                  zsigm12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
498                  zsigm22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
499                  zsigm21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
500
501                  iim1 = ji
502                  ijm1 = jj
503                  iip1 = ji + 1
504                  ijp1 = jj + 1
505                  ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   &
506                     &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) + v_ice(iim1,ijp1) )
507                  ze12 = - akappa(iim1,ijm1,2,2) * ( u_ice(iip1,ijm1) - u_ice(iim1,ijp1) - u_ice(iip1,ijp1) )   &
508                     &   - akappa(iim1,ijm1,2,1) * ( v_ice(iip1,ijm1) + v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )
509                  ze22 = - akappa(iim1,ijm1,2,2) * ( v_ice(iip1,ijm1) - v_ice(iim1,ijp1) - v_ice(iip1,ijp1) )   &
510                     &   + akappa(iim1,ijm1,2,1) * ( u_ice(iip1,ijm1) + u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )
511                  ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   &
512                     &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 
513                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
514                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
515                  zvis12 =       zviseta (iim1,ijm1) + dm
516                  zvis21 =       zviseta (iim1,ijm1)
517
518                  zdiag = zvis22 * ( ze11 + ze22 )
519                  zsigm11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
520                  zsigm12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
521                  zsigm22(ji,jj,2,2) =  zvis11 * ze22 + zdiag
522                  zsigm21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
523
524               END DO
525            END DO
526
527            ! Terms involving already up-dated velocities.     
528            !-Using the arrays zu_ice and zv_ice in the computation of the terms ze leads to JACOBI's method;
529            ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method.   
530             
531            DO jj = i_j2, i_jpjm1
532               DO ji = 2, jpim1
533                  iim1 = ji - 1
534                  ijm1 = jj - 1
535                  iip1 = ji
536                  ijp1 = jj
537                  ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) - zu_ice(iim1,ijp1) )   &
538                     &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) + zv_ice(iim1,ijp1) )
539                  ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) - zu_ice(iim1,ijp1) )   &
540                     &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) + zv_ice(iim1,ijp1) )
541                  ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) - zv_ice(iim1,ijp1) )   &
542                     &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) + zu_ice(iim1,ijp1) )
543                  ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) - zv_ice(iim1,ijp1) )   &
544                     &  -  akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) + zu_ice(iim1,ijp1) )
545                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
546                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
547                  zvis12 =       zviseta (iim1,ijm1) + dm
548                  zvis21 =       zviseta (iim1,ijm1)
549
550                  zdiag = zvis22 * ( ze11 + ze22 )
551                  zsigm11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
552                  zsigm12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
553                  zsigm22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
554                  zsigm21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
555
556
557                  iim1 = ji
558                  ijm1 = jj - 1
559                  iip1 = ji + 1
560                  ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) )   &
561                     &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) )
562                  ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )   &
563                     &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )
564                  ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )   &
565                     &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )
566                  ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) )   &
567                     &   - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) )
568                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
569                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
570                  zvis12 =       zviseta (iim1,ijm1) + dm
571                  zvis21 =       zviseta (iim1,ijm1)
572
573                  zdiag = zvis22 * ( ze11 + ze22 )
574                  zsigm11(ji,jj,2,1) =  zsigm11(ji,jj,2,1) + zvis11 * ze11 + zdiag
575                  zsigm12(ji,jj,2,1) =  zsigm12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21
576                  zsigm22(ji,jj,2,1) =  zsigm22(ji,jj,2,1) + zvis11 * ze22 + zdiag
577                  zsigm21(ji,jj,2,1) =  zsigm21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12
578
579
580                  iim1 = ji - 1
581                  ijm1 = jj 
582                  ze11 = - akappa(iim1,ijm1,1,1) * zu_ice(iim1,ijm1) + akappa(iim1,ijm1,1,2) * zv_ice(iim1,ijm1)
583                  ze12 = - akappa(iim1,ijm1,2,2) * zu_ice(iim1,ijm1) - akappa(iim1,ijm1,2,1) * zv_ice(iim1,ijm1)
584                  ze22 = - akappa(iim1,ijm1,2,2) * zv_ice(iim1,ijm1) + akappa(iim1,ijm1,2,1) * zu_ice(iim1,ijm1)
585                  ze21 = - akappa(iim1,ijm1,1,1) * zv_ice(iim1,ijm1) - akappa(iim1,ijm1,1,2) * zu_ice(iim1,ijm1)
586                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
587                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
588                  zvis12 =       zviseta (iim1,ijm1) + dm
589                  zvis21 =       zviseta (iim1,ijm1)
590
591                  zdiag = zvis22 * ( ze11 + ze22 )
592                  zsigm11(ji,jj,1,2) =  zsigm11(ji,jj,1,2) + zvis11 * ze11 + zdiag 
593                  zsigm12(ji,jj,1,2) =  zsigm12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21
594                  zsigm22(ji,jj,1,2) =  zsigm22(ji,jj,1,2) + zvis11 * ze22 + zdiag
595                  zsigm21(ji,jj,1,2) =  zsigm21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12
596
597!i             END DO
598!i          END DO
599
600!i          DO jj = i_j2, i_jpjm1
601!i             DO ji = 2, jpim1
602                  zd1(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm11(ji,jj,2,2)  &
603                     &       - alambd(ji,jj,2,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm11(ji,jj,1,2)  &
604                     &       - alambd(ji,jj,1,1,2,1) * zsigm12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm12(ji,jj,1,1)  &
605                     &       + alambd(ji,jj,1,1,2,2) * zsigm12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm12(ji,jj,1,2)  &
606                     &       + alambd(ji,jj,1,2,1,1) * zsigm21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zsigm21(ji,jj,2,1)  &
607                     &       + alambd(ji,jj,1,2,1,2) * zsigm21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zsigm21(ji,jj,2,2)  &
608                     &       - alambd(ji,jj,2,1,1,1) * zsigm22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zsigm22(ji,jj,2,1)  &
609                     &       - alambd(ji,jj,2,1,1,2) * zsigm22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zsigm22(ji,jj,2,2)
610
611                  zd2(ji,jj) = alambd(ji,jj,2,2,2,1) * zsigm21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zsigm21(ji,jj,2,2)  &
612                     &       - alambd(ji,jj,2,2,1,1) * zsigm21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zsigm21(ji,jj,1,2)  &
613                     &       - alambd(ji,jj,1,1,2,1) * zsigm22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zsigm22(ji,jj,1,1)  &
614                     &       + alambd(ji,jj,1,1,2,2) * zsigm22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zsigm22(ji,jj,1,2)  &
615                     &       - alambd(ji,jj,1,2,1,1) * zsigm11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zsigm11(ji,jj,2,1)  &
616                     &       - alambd(ji,jj,1,2,1,2) * zsigm11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zsigm11(ji,jj,2,2)  &
617                     &       + alambd(ji,jj,2,1,1,1) * zsigm12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zsigm12(ji,jj,2,1)  &
618                     &       + alambd(ji,jj,2,1,1,2) * zsigm12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zsigm12(ji,jj,2,2)
619               END DO
620            END DO
621
622            DO jj = i_j2, i_jpjm1
623               DO ji = 2, jpim1
624                  zunw = (  ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj)        &
625                     &    + ( za2(ji,jj) + zd2(ji,jj) ) * zc1(ji,jj) ) * zden(ji,jj)
626
627                  zvnw = (  ( za2(ji,jj) + zd2(ji,jj) ) * zb1(ji,jj)        &
628                     &    - ( za1(ji,jj) + zd1(ji,jj) ) * zb2(ji,jj) ) * zden(ji,jj)
629
630                  zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj)
631
632                  u_ice(ji,jj) = ( u_ice(ji,jj) + om * ( zunw - u_ice(ji,jj) ) * tmu(ji,jj) ) * zmask
633                  v_ice(ji,jj) = ( v_ice(ji,jj) + om * ( zvnw - v_ice(ji,jj) ) * tmu(ji,jj) ) * zmask
634               END DO
635            END DO
636
637            CALL lbc_lnk( u_ice, 'I', -1. )
638            CALL lbc_lnk( v_ice, 'I', -1. )
639
640            !---  5.2.5.4. Convergence test.
641            DO jj = i_j2 , i_jpjm1
642               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
643            END DO
644            zresm = MAXVAL( zresr( 1:jpi , i_j2:i_jpjm1 ) )
645
646            IF ( zresm <= resl) EXIT iflag
647
648         END DO iflag
649
650         zindu1 = 1.0 - zindu
651         DO jj = i_j1 , i_jpjm1
652            zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj)
653            zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj)
654         END DO
655      !                                                   ! ==================== !
656      END DO                                              !  end loop over iter  !
657      !                                                   ! ==================== !
658
659      IF( l_ctl .AND. lwp ) THEN
660         WRITE(numout,*) ' lim_rhg  : res= ', zresm, 'iter= ', jter,' u_ice ', SUM( u_ice ) , ' v_ice ', SUM( v_ice )
661      ENDIF
662
663   END SUBROUTINE lim_rhg
664#else
665   !!==============================================================================
666   !!                       ***  MODULE limrhg   ***
667   !!                              No sea ice
668   !!==============================================================================
669CONTAINS
670   SUBROUTINE lim_rhg         ! Empty routine
671   END SUBROUTINE lim_rhg
672
673#endif
674
675END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.