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.
limrhg.F90 in tags/start/NEMO/LIM_SRC – NEMO

source: tags/start/NEMO/LIM_SRC/limrhg.F90 @ 3489

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

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