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_2.F90 in branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 35.7 KB
Line 
1MODULE limrhg_2
2   !!======================================================================
3   !!                     ***  MODULE  limrhg_2  ***
4   !!   Ice rheology :  performs sea ice rheology
5   !!======================================================================
6   !! History :  0.0  !  1993-12  (M.A. Morales Maqueda.)  Original code
7   !!            1.0  !  1994-12  (H. Goosse)
8   !!            2.0  !  2003-12  (C. Ethe, G. Madec)  F90, mpp
9   !!             -   !  2006-08  (G. Madec)  surface module, ice-stress at I-point
10   !!             -   !  2009-09  (G. Madec)  Huge verctor optimisation
11   !!            3.3  !  2009-05  (G.Garric, C. Bricaud) addition of the lim2_evp case
12   !!----------------------------------------------------------------------
13#if defined   key_lim2   &&   defined key_lim2_vp
14   !!----------------------------------------------------------------------
15   !!   'key_lim2'                AND                   LIM-2 sea-ice model
16   !!   'key_lim2_vp'                                       VP ice rheology
17   !!----------------------------------------------------------------------
18   !!   lim_rhg_2   : computes ice velocities
19   !!----------------------------------------------------------------------
20   USE par_oce        ! ocean parameter
21   USE dom_oce        ! ocean space and time domain
22   USE sbc_oce        ! surface boundary condition: ocean variables
23   USE sbc_ice        ! surface boundary condition: ice variables
24   USE dom_ice_2      ! LIM2: ice domain
25   USE phycst         ! physical constant
26   USE ice_2          ! LIM2: ice variables
27   USE lbclnk         ! lateral boundary condition - MPP exchanges
28   USE lib_mpp        ! MPP library
29   USE wrk_nemo       ! work arrays
30   USE in_out_manager ! I/O manager
31   USE prtctl         ! Print control
32   USE oce     , ONLY : snwice_mass, snwice_mass_b
33   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
34#if defined key_agrif
35   USE agrif_lim2_interp ! nesting
36#endif
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   lim_rhg_2         ! routine called by lim_dyn
42
43   REAL(wp) ::   rzero   = 0._wp   ! constant value: zero
44   REAL(wp) ::   rone    = 1._wp   !            and  one
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE lim_rhg_2( k_j1, k_jpj )
56      !!-------------------------------------------------------------------
57      !!                 ***  SUBROUTINR lim_rhg_2  ***
58      !!
59      !! ** purpose :   determines the velocity field of sea ice by using
60      !!  atmospheric (wind stress) and oceanic (water stress and surface
61      !!  tilt) forcings. Ice-ice interaction is described by a non-linear
62      !!  viscous-plastic law including shear strength and a bulk rheology.
63      !!
64      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity defined
65      !!              at I-point
66      !!-------------------------------------------------------------------
67      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
68      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
69      !!
70      INTEGER ::   ji, jj              ! dummy loop indices
71      INTEGER ::   iter, jter          ! temporary integers
72      CHARACTER (len=50) ::   charout
73      REAL(wp) ::   ze11  , ze12  , ze22  , ze21               ! local scalars
74      REAL(wp) ::   zt11  , zt12  , zt21  , zt22               !   -      -
75      REAL(wp) ::   zvis11, zvis21, zvis12, zvis22             !   -      -
76      REAL(wp) ::   zgphsx, ztagnx, zgsshx, zunw, zur, zusw    !   -      -
77      REAL(wp) ::   zgphsy, ztagny, zgsshy, zvnw, zvr          !   -      -
78      REAL(wp) ::   zresm,  za, zac, zmod
79      REAL(wp) ::   zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal
80      REAL(wp) ::   ztrace2, zdeter, zdelta, zmask, zdgp, zdgi, zdiag
81      REAL(wp) ::   za1, zb1, zc1, zd1
82      REAL(wp) ::   za2, zb2, zc2, zd2, zden
83      REAL(wp) ::   zs11_11, zs11_12, zs11_21, zs11_22
84      REAL(wp) ::   zs12_11, zs12_12, zs12_21, zs12_22
85      REAL(wp) ::   zs21_11, zs21_12, zs21_21, zs21_22
86      REAL(wp) ::   zs22_11, zs22_12, zs22_21, zs22_22
87      REAL(wp) ::   zintb, zintn
88      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld, zmass, zcorl
89      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct, za2ct, zresr
90      REAL(wp), POINTER, DIMENSION(:,:) ::   zc1u, zc1v, zc2u, zc2v
91      REAL(wp), POINTER, DIMENSION(:,:) ::   zsang, zpice
92      REAL(wp), POINTER, DIMENSION(:,:) ::   zu0, zv0
93      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_n, zv_n
94      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_a, zv_a
95      REAL(wp), POINTER, DIMENSION(:,:) ::   zviszeta, zviseta
96      REAL(wp), POINTER, DIMENSION(:,:) ::   zzfrld, zztms
97      REAL(wp), POINTER, DIMENSION(:,:) ::   zi1, zi2, zmasst, zpresh
98      !!-------------------------------------------------------------------
99     
100      CALL wrk_alloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr )
101      CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice )
102      CALL wrk_alloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 )
103      CALL wrk_alloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 )
104
105      !  Store initial velocities
106      !  ----------------
107      zztms(:,0    ) = 0._wp        ;   zzfrld(:,0    ) = 0._wp
108      zztms(:,jpj+1) = 0._wp        ;   zzfrld(:,jpj+1) = 0._wp
109      zu0  (:,0    ) = 0._wp        ;   zv0   (:,0    ) = 0._wp
110      zu0  (:,jpj+1) = 0._wp        ;   zv0   (:,jpj+1) = 0._wp
111      zztms(:,1:jpj) = tms  (:,:)   ;   zzfrld(:,1:jpj) = frld (:,:)
112      zu0  (:,1:jpj) = u_ice(:,:)   ;   zv0   (:,1:jpj) = v_ice(:,:)
113      zu_a (:, :   ) = zu0  (:,:)   ;   zv_a  (:, :   ) = zv0  (:,:)
114      zu_n (:, :   ) = zu0  (:,:)   ;   zv_n  (:, :   ) = zv0  (:,:)
115
116!i
117      zi1   (:,:) = 0._wp
118      zi2   (:,:) = 0._wp
119      zpresh(:,:) = 0._wp
120      zmasst(:,:) = 0._wp
121!i
122!!gm violant
123      zfrld(:,:) =0._wp
124      zcorl(:,:) =0._wp
125      zmass(:,:) =0._wp
126      za1ct(:,:) =0._wp
127      za2ct(:,:) =0._wp
128!!gm end
129
130      zviszeta(:,:) = 0._wp
131      zviseta (:,:) = 0._wp
132
133!i    zviszeta(:,0    ) = 0._wp    ;    zviseta(:,0    ) = 0._wp
134!i    zviszeta(:,jpj  ) = 0._wp    ;    zviseta(:,jpj  ) = 0._wp
135!i    zviszeta(:,jpj+1) = 0._wp    ;    zviseta(:,jpj+1) = 0._wp
136
137      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
138          !
139          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
140          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
141         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
142          !
143          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
144          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
145         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
146          !
147         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
148          !
149         !
150      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
151         zpice(:,:) = ssh_m(:,:)
152      ENDIF
153#if defined key_agrif
154      ! load the boundary value of velocity in special array zuive and zvice
155      CALL agrif_rhg_lim2_load
156#endif
157
158      ! Ice mass, ice strength, and wind stress at the center            |
159      ! of the grid squares.                                             |
160      !-------------------------------------------------------------------
161
162      DO jj = k_j1 , k_jpj-1
163         DO ji = 1 , jpi
164            ! only the sinus changes its sign with the hemisphere
165            zsang(ji,jj)  = SIGN( 1._wp, fcor(ji,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere
166            !
167            zmasst(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
168            zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
169!!gm  :: stress given at I-point (F-point for the ocean) only compute the ponderation with the ice fraction (1-frld)
170            zi1(ji,jj)    = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
171            zi2(ji,jj)    = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
172         END DO
173      END DO
174
175
176      !---------------------------------------------------------------------------
177      !  Wind stress, coriolis and mass terms at the corners of the grid squares |
178      !  Gradient of ice strenght.                                               |
179      !---------------------------------------------------------------------------
180         
181      DO jj = k_j1+1, k_jpj-1
182         DO ji = 2, jpi    ! NO vector opt.
183            zstms = zztms(ji,jj  ) * wght(ji,jj,2,2) + zztms(ji-1,jj  ) * wght(ji,jj,1,2)   &
184               &  + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1)
185            zusw  = 1._wp / MAX( zstms, epsd )
186
187            zt11 = zztms(ji  ,jj  ) * zzfrld(ji  ,jj  ) 
188            zt12 = zztms(ji-1,jj  ) * zzfrld(ji-1,jj  ) 
189            zt21 = zztms(ji  ,jj-1) * zzfrld(ji  ,jj-1) 
190            zt22 = zztms(ji-1,jj-1) * zzfrld(ji-1,jj-1)
191
192            ! Leads area.
193            zfrld(ji,jj) =  (  zt11 * wght(ji,jj,2,2) + zt12 * wght(ji,jj,1,2)   &
194               &             + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw
195
196            ! Mass and coriolis coeff. at I-point
197            zmass(ji,jj) = ( zmasst(ji,jj  ) * wght(ji,jj,2,2) + zmasst(ji-1,jj  ) * wght(ji,jj,1,2)   &
198               &           + zmasst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
199            zcorl(ji,jj) = zmass(ji,jj) &
200               &           *( fcor(ji,jj  ) * wght(ji,jj,2,2) + fcor(ji-1,jj  )*wght(ji,jj,1,2)   &
201               &            + fcor(ji,jj-1) * wght(ji,jj,2,1) + fcor(ji-1,jj-1)*wght(ji,jj,1,1) ) * zusw
202
203            ! Wind stress.
204            ! always provide stress at I-point
205            ztagnx = ( zi1(ji,jj  ) * wght(ji,jj,2,2) + zi1(ji-1,jj  ) * wght(ji,jj,1,2)   &
206               &     + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utau_ice(ji,jj)
207            ztagny = ( zi2(ji,jj  ) * wght(ji,jj,2,2) + zi2(ji-1,jj  ) * wght(ji,jj,1,2)   &
208               &     + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtau_ice(ji,jj)
209
210            ! Gradient of ice strength
211            zgphsx =   ( alambd(ji,jj,2,2,2,1) - alambd(ji,jj,2,1,2,1) ) * zpresh(ji  ,jj-1)   &
212               &     + ( alambd(ji,jj,2,2,2,2) - alambd(ji,jj,2,1,2,2) ) * zpresh(ji  ,jj  )   &
213               &     - ( alambd(ji,jj,2,2,1,1) + alambd(ji,jj,2,1,1,1) ) * zpresh(ji-1,jj-1)   &
214               &     - ( alambd(ji,jj,2,2,1,2) + alambd(ji,jj,2,1,1,2) ) * zpresh(ji-1,jj  )
215
216            zgphsy = - ( alambd(ji,jj,1,1,2,1) + alambd(ji,jj,1,2,2,1) ) * zpresh(ji  ,jj-1)   &
217               &     - ( alambd(ji,jj,1,1,1,1) + alambd(ji,jj,1,2,1,1) ) * zpresh(ji-1,jj-1)   &
218               &     + ( alambd(ji,jj,1,1,2,2) - alambd(ji,jj,1,2,2,2) ) * zpresh(ji  ,jj  )   &
219               &     + ( alambd(ji,jj,1,1,1,2) - alambd(ji,jj,1,2,1,2) ) * zpresh(ji-1,jj  )
220
221            ! Gradient of the sea surface height
222            zgsshx =  (   (zpice(ji  ,jj  ) - zpice(ji-1,jj  ))/e1u(ji-1,jj  )   &
223               &       +  (zpice(ji  ,jj-1) - zpice(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp
224            zgsshy =  (   (zpice(ji  ,jj  ) - zpice(ji  ,jj-1))/e2v(ji  ,jj-1)   &
225               &       +  (zpice(ji-1,jj  ) - zpice(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp
226
227            ! Computation of the velocity field taking into account the ice-ice interaction.                                 
228            ! Terms that are independent of the ice velocity field.
229            za1ct(ji,jj) = ztagnx - zmass(ji,jj) * grav * zgsshx - zgphsx
230            za2ct(ji,jj) = ztagny - zmass(ji,jj) * grav * zgsshy - zgphsy
231         END DO
232      END DO
233
234
235      ! SOLUTION OF THE MOMENTUM EQUATION.
236      !------------------------------------------
237      !                                                   ! ==================== !
238      DO iter = 1 , 2 * nbiter                            !    loop over iter    !
239         !                                                ! ==================== !       
240         zindu = MOD( iter , 2 )
241         zusdtp = ( zindu * 2._wp + ( 1._wp - zindu ) * 1._wp )  * REAL( nbiter ) / rdt_ice
242
243         ! Computation of free drift field for free slip boundary conditions.
244
245         DO jj = k_j1, k_jpj-1
246            DO ji = 1, fs_jpim1
247               !- Rate of strain tensor.
248               zt11 =   akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji,jj  ) - zu_a(ji  ,jj+1) )  &
249                  &   + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji,jj  ) + zv_a(ji  ,jj+1) )
250               zt12 = - akappa(ji,jj,2,2) * ( zu_a(ji  ,jj) + zu_a(ji+1,jj  ) - zu_a(ji,jj+1) - zu_a(ji+1,jj+1) )  &
251                  &   - akappa(ji,jj,2,1) * ( zv_a(ji  ,jj) + zv_a(ji+1,jj  ) + zv_a(ji,jj+1) + zv_a(ji+1,jj+1) )
252               zt22 = - akappa(ji,jj,2,2) * ( zv_a(ji  ,jj) + zv_a(ji+1,jj  ) - zv_a(ji,jj+1) - zv_a(ji+1,jj+1) )  &
253                  &   + akappa(ji,jj,2,1) * ( zu_a(ji  ,jj) + zu_a(ji+1,jj  ) + zu_a(ji,jj+1) + zu_a(ji+1,jj+1) )
254               zt21 =   akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji,jj  ) - zv_a(ji  ,jj+1) )  &
255                  &   - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji,jj  ) + zu_a(ji  ,jj+1) )
256
257               !- Rate of strain tensor.
258               zdgp = zt11 + zt22
259               zdgi = zt12 + zt21
260               ztrace2 = zdgp * zdgp 
261               zdeter  = zt11 * zt22 - 0.25_wp * zdgi * zdgi
262
263               !  Creep limit depends on the size of the grid.
264               zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4._wp * zdeter ) * usecc2 ),  rn_creepl)
265
266               !-  Computation of viscosities.
267               zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn )
268               zviseta (ji,jj) = zviszeta(ji,jj) * usecc2
269            END DO
270         END DO
271
272         !-  Determination of zc1u, zc2u, zc1v and zc2v.
273         DO jj = k_j1+1, k_jpj-1
274            DO ji = 2, fs_jpim1   ! NO vector opt.
275               !* zc1u , zc2v
276               zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
277               zvis12 =         zviseta (ji-1,jj-1) + dm
278               zvis21 =         zviseta (ji-1,jj-1)
279               zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
280               zdiag  = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) )
281               zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,1) + zdiag
282               zs12_11 =  zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2)
283               zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2)
284               zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,1) + zdiag
285
286               zvis11 = 2._wp * zviseta (ji,jj-1) + dm
287               zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
288               zvis12 =         zviseta (ji,jj-1) + dm
289               zvis21 =         zviseta (ji,jj-1)
290               zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) )
291               zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag
292               zs12_21 =  zvis12 * akappa(ji,jj-1,2,2) - zvis21 * akappa(ji,jj-1,1,2)
293               zs22_21 =  zvis11 * akappa(ji,jj-1,2,1) + zdiag
294               zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2)
295
296               zvis11 = 2._wp * zviseta (ji-1,jj) + dm
297               zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
298               zvis12 =         zviseta (ji-1,jj) + dm
299               zvis21 =         zviseta (ji-1,jj)
300               zdiag  = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) )
301               zs11_12 =  zvis11 * akappa(ji-1,jj,1,1) + zdiag
302               zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2)
303               zs22_12 =  zvis11 * akappa(ji-1,jj,2,1) + zdiag
304               zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2)
305
306               zvis11 = 2._wp * zviseta (ji,jj) + dm
307               zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
308               zvis12 =         zviseta (ji,jj) + dm
309               zvis21 =         zviseta (ji,jj)
310               zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) )
311               zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag
312               zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2)
313               zs22_22 =  zvis11 * akappa(ji,jj,2,1) + zdiag
314               zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2)
315
316               zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22   &
317                  &          - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12   &
318                  &          - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11   &
319                  &          + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12   &
320                  &          + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21   &
321                  &          + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22   &
322                  &          - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21   &
323                  &          - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
324
325               zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22   &
326                  &          - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12   &
327                  &          - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11   &
328                  &          + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12   &
329                  &          - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21   &
330                  &          - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22   &
331                  &          + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21   &
332                  &          + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
333
334               !* zc1v , zc2v.
335               zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
336               zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
337               zvis12 =         zviseta (ji-1,jj-1) + dm
338               zvis21 =         zviseta (ji-1,jj-1)
339               zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) )
340               zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,2) + zdiag
341               zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1)
342               zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,2) + zdiag
343               zs21_11 =  zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1)
344 
345               zvis11 = 2._wp * zviseta (ji,jj-1) + dm
346               zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
347               zvis12 =         zviseta (ji,jj-1) + dm
348               zvis21 =         zviseta (ji,jj-1)
349               zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) )
350               zs11_21 =  zvis11 * akappa(ji,jj-1,1,2) + zdiag
351               zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1)
352               zs22_21 =  zvis11 * akappa(ji,jj-1,2,2) + zdiag
353               zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1)
354
355               zvis11 = 2._wp * zviseta (ji-1,jj) + dm
356               zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
357               zvis12 =         zviseta (ji-1,jj) + dm
358               zvis21 =         zviseta (ji-1,jj)
359               zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) )
360               zs11_12 =  zvis11 * akappa(ji-1,jj,1,2) + zdiag
361               zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1)
362               zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag
363               zs21_12 =  zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1)
364
365               zvis11 = 2._wp * zviseta (ji,jj) + dm
366               zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
367               zvis12 =         zviseta (ji,jj) + dm
368               zvis21 =         zviseta (ji,jj)
369               zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) )
370               zs11_22 =  zvis11 * akappa(ji,jj,1,2) + zdiag
371               zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1)
372               zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag
373               zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1)
374
375               zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22   &
376                  &          - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12   &
377                  &          - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11   &
378                  &          + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12   &
379                  &          + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21   &
380                  &          + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22   &
381                  &          - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21   &
382                  &          - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
383
384               zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22   &
385                  &          - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12   &
386                  &          - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11   &
387                  &          + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12   &
388                  &          - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21   &
389                  &          - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22   &
390                  &          + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21   &
391                  &          + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
392            END DO
393         END DO
394
395         ! GAUSS-SEIDEL method
396         !                                                      ! ================ !
397iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    !
398            !                                                   ! ================ !
399            DO jj = k_j1+1, k_jpj-1
400               DO ji = 2, fs_jpim1   ! NO vector opt.
401                  !
402                  ze11 =   akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj)
403                  ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj)
404                  ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj)
405                  ze21 =   akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj)
406                  zvis11 = 2._wp * zviseta (ji,jj-1) + dm
407                  zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
408                  zvis12 =         zviseta (ji,jj-1) + dm
409                  zvis21 =         zviseta (ji,jj-1)
410                  zdiag = zvis22 * ( ze11 + ze22 )
411                  zs11_21 =  zvis11 * ze11 + zdiag
412                  zs12_21 =  zvis12 * ze12 + zvis21 * ze21
413                  zs22_21 =  zvis11 * ze22 + zdiag
414                  zs21_21 =  zvis12 * ze21 + zvis21 * ze12
415
416                  ze11 =   akappa(ji-1,jj,1,1) * ( zu_a(ji  ,jj+1) - zu_a(ji-1,jj+1) )   &
417                     &   + akappa(ji-1,jj,1,2) * ( zv_a(ji  ,jj+1) + zv_a(ji-1,jj+1) )
418                  ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji  ,jj+1) )   &
419                     &   - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji  ,jj+1) )
420                  ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji  ,jj+1) )   &
421                     &   + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji  ,jj+1) )
422                  ze21 =   akappa(ji-1,jj,1,1) * ( zv_a(ji  ,jj+1) - zv_a(ji-1,jj+1) )   &
423                     &   - akappa(ji-1,jj,1,2) * ( zu_a(ji  ,jj+1) + zu_a(ji-1,jj+1) )
424                  zvis11 = 2._wp * zviseta (ji-1,jj) + dm
425                  zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
426                  zvis12 =         zviseta (ji-1,jj) + dm
427                  zvis21 =         zviseta (ji-1,jj)
428                  zdiag = zvis22 * ( ze11 + ze22 )
429                  zs11_12 =  zvis11 * ze11 + zdiag
430                  zs12_12 =  zvis12 * ze12 + zvis21 * ze21
431                  zs22_12 =  zvis11 * ze22 + zdiag
432                  zs21_12 =  zvis12 * ze21 + zvis21 * ze12
433
434                  ze11 =   akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji  ,jj+1) )   &
435                     &   + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji  ,jj+1) )
436                  ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji  ,jj+1) - zu_a(ji+1,jj+1) )   &
437                     &   - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji  ,jj+1) + zv_a(ji+1,jj+1) )
438                  ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji  ,jj+1) - zv_a(ji+1,jj+1) )   &
439                     &   + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji  ,jj+1) + zu_a(ji+1,jj+1) )
440                  ze21 =   akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji  ,jj+1) )   &
441                     &   - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji  ,jj+1) )
442                  zvis11 = 2._wp * zviseta (ji,jj) + dm
443                  zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
444                  zvis12 =         zviseta (ji,jj) + dm
445                  zvis21 =         zviseta (ji,jj)
446                  zdiag = zvis22 * ( ze11 + ze22 )
447                  zs11_22 =  zvis11 * ze11 + zdiag
448                  zs12_22 =  zvis12 * ze12 + zvis21 * ze21
449                  zs22_22 =  zvis11 * ze22 + zdiag
450                  zs21_22 =  zvis12 * ze21 + zvis21 * ze12
451
452            ! 2nd part
453                  ze11 =   akappa(ji-1,jj-1,1,1) * ( zu_a(ji  ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) )   &
454                     &   + akappa(ji-1,jj-1,1,2) * ( zv_a(ji  ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) )
455                  ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji  ,jj-1) - zu_a(ji-1,jj) )   &
456                     &   - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji  ,jj-1) + zv_a(ji-1,jj) )
457                  ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji  ,jj-1) - zv_a(ji-1,jj) )   &
458                     &   + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji  ,jj-1) + zu_a(ji-1,jj) )
459                  ze21 =   akappa(ji-1,jj-1,1,1) * ( zv_a(ji  ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) )   &
460                     &  -  akappa(ji-1,jj-1,1,2) * ( zu_a(ji  ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) )
461                  zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
462                  zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
463                  zvis12 =         zviseta (ji-1,jj-1) + dm
464                  zvis21 =         zviseta (ji-1,jj-1)
465                  zdiag = zvis22 * ( ze11 + ze22 )
466                  zs11_11 =  zvis11 * ze11 + zdiag
467                  zs12_11 =  zvis12 * ze12 + zvis21 * ze21
468                  zs22_11 =  zvis11 * ze22 + zdiag
469                  zs21_11 =  zvis12 * ze21 + zvis21 * ze12
470
471                  ze11 =   akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji  ,jj-1) )   &
472                     &   + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji  ,jj-1) )
473                  ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji  ,jj-1) + zu_a(ji+1,jj-1) )   &
474                     &   - akappa(ji,jj-1,2,1) * ( zv_a(ji  ,jj-1) + zv_a(ji+1,jj-1) )
475                  ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji  ,jj-1) + zv_a(ji+1,jj-1) )   &
476                     &   + akappa(ji,jj-1,2,1) * ( zu_a(ji  ,jj-1) + zu_a(ji+1,jj-1) )
477                  ze21 =   akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji  ,jj-1) )   &
478                     &   - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji  ,jj-1) )
479                  zvis11 = 2._wp * zviseta (ji,jj-1) + dm
480                  zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
481                  zvis12 =         zviseta (ji,jj-1) + dm
482                  zvis21 =         zviseta (ji,jj-1)
483                  zdiag = zvis22 * ( ze11 + ze22 )
484                  zs11_21 =  zs11_21 + zvis11 * ze11 + zdiag
485                  zs12_21 =  zs12_21 + zvis12 * ze12 + zvis21 * ze21
486                  zs22_21 =  zs22_21 + zvis11 * ze22 + zdiag
487                  zs21_21 =  zs21_21 + zvis12 * ze21 + zvis21 * ze12
488
489                  ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj)
490                  ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj)
491                  ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj)
492                  ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj)
493                  zvis11 = 2._wp * zviseta (ji-1,jj) + dm
494                  zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
495                  zvis12 =         zviseta (ji-1,jj) + dm
496                  zvis21 =         zviseta (ji-1,jj)
497                  zdiag = zvis22 * ( ze11 + ze22 )
498                  zs11_12 =  zs11_12 + zvis11 * ze11 + zdiag
499                  zs12_12 =  zs12_12 + zvis12 * ze12 + zvis21 * ze21
500                  zs22_12 =  zs22_12 + zvis11 * ze22 + zdiag
501                  zs21_12 =  zs21_12 + zvis12 * ze21 + zvis21 * ze12
502
503                  zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22  &
504                     &  - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12  &
505                     &  - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11  &
506                     &  + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12  &
507                     &  + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21  &
508                     &  + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22  &
509                     &  - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21  &
510                     &  - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
511
512                  zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22  &
513                     &  - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12  &
514                     &  - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11  &
515                     &  + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12  &
516                     &  - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21  &
517                     &  - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22  &
518                     &  + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21  &
519                     &  + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
520
521                  zur     = zu_a(ji,jj) - u_oce(ji,jj)
522                  zvr     = zv_a(ji,jj) - v_oce(ji,jj)
523!!!!
524                  zmod    = SQRT( zur*zur + zvr*zvr ) * ( 1._wp - zfrld(ji,jj) )
525                  za      = rhoco * zmod
526!!!!
527!!gm chg resul    za      = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1._wp - zfrld(ji,jj) )
528                  zac     = za * cangvg
529                  zmpzas  = alpha * zcorl(ji,jj) + za * zsang(ji,jj)
530                  zmassdt = zusdtp * zmass(ji,jj)
531                  zcorlal = ( 1._wp - alpha ) * zcorl(ji,jj)
532
533                  za1 =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   &
534                     &        + za * ( cangvg * u_oce(ji,jj) - zsang(ji,jj) * v_oce(ji,jj) )
535                  za2 =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   &
536                     &        + za * ( cangvg * v_oce(ji,jj) + zsang(ji,jj) * u_oce(ji,jj) )
537                  zb1    = zmassdt + zac - zc1u(ji,jj)
538                  zb2    = zmpzas        - zc2u(ji,jj)
539                  zc1    = zmpzas        + zc1v(ji,jj)
540                  zc2    = zmassdt + zac - zc2v(ji,jj)
541                  zdeter = zc1 * zb2 + zc2 * zb1
542                  zden   = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) )
543                  zunw   = (  ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden
544                  zvnw   = (  ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden
545                  zmask  = ( 1._wp - MAX( rzero, SIGN( rone , 1._wp - zmass(ji,jj) ) ) ) * tmu(ji,jj)
546
547                  zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask
548                  zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask
549               END DO
550            END DO
551
552            CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. )
553            CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. )
554
555#if defined key_agrif
556            ! copy the boundary value from u_ice_nst and v_ice_nst to u_ice and v_ice
557            ! before next interations
558            CALL agrif_rhg_lim2(zu_n,zv_n)
559#endif
560
561            ! Test of Convergence
562            DO jj = k_j1+1 , k_jpj-1
563               zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) )
564            END DO
565            zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) )
566!!!! this should be faster on scalar processor
567!           zresm = MAXVAL(  MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ),   &
568!              &                  ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) )  )
569!!!!
570            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
571
572            DO jj = k_j1, k_jpj
573               zu_a(:,jj) = zu_n(:,jj)
574               zv_a(:,jj) = zv_n(:,jj)
575            END DO
576
577            IF( zresm <= resl )   EXIT   iflag
578
579            !                                                   ! ================ !
580         END DO    iflag                                        !  end Relaxation  !
581         !                                                      ! ================ !
582
583         IF( zindu == 0 ) THEN      ! even iteration
584            DO jj = k_j1 , k_jpj-1
585               zu0(:,jj) = zu_a(:,jj)
586               zv0(:,jj) = zv_a(:,jj)
587            END DO
588         ENDIF
589         !                                                ! ==================== !
590      END DO                                              !  end loop over iter  !
591      !                                                   ! ==================== !
592
593      u_ice(:,:) = zu_a(:,1:jpj)
594      v_ice(:,:) = zv_a(:,1:jpj)
595
596      IF(ln_ctl) THEN
597         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
598         CALL prt_ctl_info(charout)
599         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
600      ENDIF
601
602      CALL wrk_dealloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr )
603      CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice )
604      CALL wrk_dealloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 )
605      CALL wrk_dealloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 )
606
607   END SUBROUTINE lim_rhg_2
608
609#else
610   !!----------------------------------------------------------------------
611   !!   Default option        Dummy module      NO VP & LIM-2 sea-ice model
612   !!----------------------------------------------------------------------
613CONTAINS
614   SUBROUTINE lim_rhg_2( k1 , k2 )       ! Dummy routine
615      WRITE(*,*) 'lim_rhg_2: You should not have seen this print! error?', k1, k2
616   END SUBROUTINE lim_rhg_2
617#endif
618
619   !!==============================================================================
620END MODULE limrhg_2
Note: See TracBrowser for help on using the repository browser.