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/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

source: branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90 @ 6808

Last change on this file since 6808 was 6808, checked in by jamesharle, 8 years ago

merge with trunk@6232 for consistency with SSB code

  • 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.