source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 7753

Last change on this file since 7753 was 7753, checked in by mocavero, 4 years ago

Reverting trunk to remove OpenMP

  • Property svn:keywords set to Id
File size: 39.2 KB
Line 
1MODULE limrhg
2   !!======================================================================
3   !!                     ***  MODULE  limrhg  ***
4   !!   Ice rheology : sea ice rheology
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
8   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
9   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas
10   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation
11   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
12   !!            3.6  !  2016-06  (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
13   !!----------------------------------------------------------------------
14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                      LIM-3 sea-ice model
17   !!----------------------------------------------------------------------
18   !!   lim_rhg       : computes ice velocities
19   !!----------------------------------------------------------------------
20   USE phycst         ! Physical constant
21   USE oce     , ONLY :  snwice_mass, snwice_mass_b
22   USE par_oce        ! Ocean parameters
23   USE dom_oce        ! Ocean domain
24   USE sbc_oce        ! Surface boundary condition: ocean fields
25   USE sbc_ice        ! Surface boundary condition: ice fields
26   USE ice            ! ice variables
27   USE limitd_me      ! ice strength
28   USE lbclnk         ! Lateral Boundary Condition / MPP link
29   USE lib_mpp        ! MPP library
30   USE wrk_nemo       ! work arrays
31   USE in_out_manager ! I/O manager
32   USE prtctl         ! Print control
33   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
34#if defined key_agrif
35   USE agrif_lim3_interp
36#endif
37   USE bdy_oce   , ONLY: ln_bdy 
38   USE bdyice_lim 
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   lim_rhg        ! routine called by lim_dyn
44
45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE lim_rhg
55      !!-------------------------------------------------------------------
56      !!                 ***  SUBROUTINE lim_rhg  ***
57      !!                          EVP-C-grid
58      !!
59      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
60      !!  stress and sea-surface slope. Ice-ice interaction is described by
61      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
62      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
63      !!
64      !!  The points in the C-grid look like this, dear reader
65      !!
66      !!                              (ji,jj)
67      !!                                 |
68      !!                                 |
69      !!                      (ji-1,jj)  |  (ji,jj)
70      !!                             ---------   
71      !!                            |         |
72      !!                            | (ji,jj) |------(ji,jj)
73      !!                            |         |
74      !!                             ---------   
75      !!                     (ji-1,jj-1)     (ji,jj-1)
76      !!
77      !! ** Inputs  : - wind forcing (stress), oceanic currents
78      !!                ice total volume (vt_i) per unit area
79      !!                snow total volume (vt_s) per unit area
80      !!
81      !! ** Action  : - compute u_ice, v_ice : the components of the
82      !!                sea-ice velocity vector
83      !!              - compute delta_i, shear_i, divu_i, which are inputs
84      !!                of the ice thickness distribution
85      !!
86      !! ** Steps   : 1) Compute ice snow mass, ice strength
87      !!              2) Compute wind, oceanic stresses, mass terms and
88      !!                 coriolis terms of the momentum equation
89      !!              3) Solve the momentum equation (iterative procedure)
90      !!              4) Prevent high velocities if the ice is thin
91      !!              5) Recompute invariants of the strain rate tensor
92      !!                 which are inputs of the ITD, store stress
93      !!                 for the next time step
94      !!              6) Control prints of residual (convergence)
95      !!                 and charge ellipse.
96      !!                 The user should make sure that the parameters
97      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
98      !!                 on the charge ellipse for plastic flow
99      !!                 e.g. in the Canadian Archipelago
100      !!
101      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
102      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
103      !!              but this solution appears very unstable (see Kimmritz et al 2016)
104      !!
105      !! References : Hunke and Dukowicz, JPO97
106      !!              Bouillon et al., Ocean Modelling 2009
107      !!              Bouillon et al., Ocean Modelling 2013
108      !!-------------------------------------------------------------------
109      INTEGER ::   ji, jj       ! dummy loop indices
110      INTEGER ::   jter         ! local integers
111      CHARACTER (len=50) ::   charout
112
113      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
114      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
115      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
116      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
117      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
118      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
119      REAL(wp) ::   zTauO, zTauB, zTauE, zCor, zvel                          ! temporary scalars
120
121      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
122      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
123      REAL(wp) ::   zintb, zintn                                             ! dummy argument
124     
125      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors
126      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                         ! P/delta at T points
127      !
128      REAL(wp), POINTER, DIMENSION(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points
129      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
130      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points
131      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
132      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
133      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
134      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses
135     
136      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! shear
137      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12                  ! stress tensor components
138      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr           ! check convergence
139      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                           ! array used for the calculation of ice surface slope:
140                                                                             !   ocean surface (ssh_m) if ice is not embedded
141                                                                             !   ice top surface if ice is embedded   
142      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays
143      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
144      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice
145
146      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
147      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
148      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
149      !!-------------------------------------------------------------------
150
151      CALL wrk_alloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
152      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
153      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
154      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
155      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
156
157#if defined key_agrif 
158      CALL agrif_interp_lim3( 'U', 0, nn_nevp )   ! First interpolation of coarse values
159      CALL agrif_interp_lim3( 'V', 0, nn_nevp )
160#endif
161      !
162      !------------------------------------------------------------------------------!
163      ! 0) mask at F points for the ice
164      !------------------------------------------------------------------------------!
165      ! ocean/land mask
166      DO jj = 1, jpjm1
167         DO ji = 1, jpim1      ! NO vector opt.
168            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
169         END DO
170      END DO
171      CALL lbc_lnk( zfmask, 'F', 1._wp )
172
173      ! Lateral boundary conditions on velocity (modify zfmask)
174      zwf(:,:) = zfmask(:,:)
175      DO jj = 2, jpjm1
176         DO ji = fs_2, fs_jpim1   ! vector opt.
177            IF( zfmask(ji,jj) == 0._wp ) THEN
178               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
179            ENDIF
180         END DO
181      END DO
182      DO jj = 2, jpjm1
183         IF( zfmask(1,jj) == 0._wp ) THEN
184            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
185         ENDIF
186         IF( zfmask(jpi,jj) == 0._wp ) THEN
187            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
188         ENDIF
189      END DO
190      DO ji = 2, jpim1
191         IF( zfmask(ji,1) == 0._wp ) THEN
192            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
193         ENDIF
194         IF( zfmask(ji,jpj) == 0._wp ) THEN
195            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
196         ENDIF
197      END DO
198      CALL lbc_lnk( zfmask, 'F', 1._wp )
199
200      !------------------------------------------------------------------------------!
201      ! 1) define some variables and initialize arrays
202      !------------------------------------------------------------------------------!
203      zrhoco = rau0 * rn_cio 
204
205      ! ecc2: square of yield ellipse eccenticrity
206      ecc2    = rn_ecc * rn_ecc
207      z1_ecc2 = 1._wp / ecc2
208
209      ! Time step for subcycling
210      zdtevp   = rdt_ice / REAL( nn_nevp )
211      z1_dtevp = 1._wp / zdtevp
212
213      ! alpha parameters (Bouillon 2009)
214      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
215      zalph2 = zalph1 * z1_ecc2
216
217      ! alpha and beta parameters (Bouillon 2013)
218      !!zalph1 = 40.
219      !!zalph2 = 40.
220      !!zbeta  = 3000.
221      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001)
222
223      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
224      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
225
226      ! Initialise stress tensor
227      zs1 (:,:) = stress1_i (:,:) 
228      zs2 (:,:) = stress2_i (:,:)
229      zs12(:,:) = stress12_i(:,:)
230
231      ! Ice strength
232      CALL lim_itd_me_icestrength( nn_icestr )
233
234      ! scale factors
235      DO jj = 2, jpjm1
236         DO ji = fs_2, fs_jpim1
237            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
238            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
239         END DO
240      END DO
241           
242      !
243      !------------------------------------------------------------------------------!
244      ! 2) Wind / ocean stress, mass terms, coriolis terms
245      !------------------------------------------------------------------------------!
246
247      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
248         !                                           
249         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
250         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
251         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
252         !
253         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
254         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
255         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
256         !
257         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
258         !
259      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
260         zpice(:,:) = ssh_m(:,:)
261      ENDIF
262
263      DO jj = 2, jpjm1
264         DO ji = fs_2, fs_jpim1
265
266            ! ice fraction at U-V points
267            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
268            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
269
270            ! Ice/snow mass at U-V points
271            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
272            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
273            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
274            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
275            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
276
277            ! Ocean currents at U-V points
278            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
279               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
280           
281            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
282               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
283
284            ! Coriolis at T points (m*f)
285            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
286
287            ! m/dt
288            zmU_t(ji,jj)    = zmassU * z1_dtevp
289            zmV_t(ji,jj)    = zmassV * z1_dtevp
290
291            ! Drag ice-atm.
292            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
293            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
294
295            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
296            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
297            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
298
299            ! masks
300            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
301            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
302
303            ! switches
304            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
305            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
306
307         END DO
308      END DO
309      CALL lbc_lnk( zmf, 'T', 1. )
310      !
311      !------------------------------------------------------------------------------!
312      ! 3) Solution of the momentum equation, iterative procedure
313      !------------------------------------------------------------------------------!
314      !
315      !                                               !----------------------!
316      DO jter = 1 , nn_nevp                           !    loop over jter    !
317         !                                            !----------------------!       
318         IF(ln_ctl) THEN   ! Convergence test
319            DO jj = 1, jpjm1
320               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
321               zv_ice(:,jj) = v_ice(:,jj)
322            END DO
323         ENDIF
324
325         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
326         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
327            DO ji = 1, jpim1
328
329               ! shear at F points
330               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
331                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
332                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
333
334            END DO
335         END DO
336         CALL lbc_lnk( zds, 'F', 1. )
337
338         DO jj = 2, jpjm1
339            DO ji = 2, jpim1 ! no vector loop
340
341               ! shear**2 at T points (doc eq. A16)
342               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
343                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
344                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
345             
346               ! divergence at T points
347               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
348                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
349                  &    ) * r1_e1e2t(ji,jj)
350               zdiv2 = zdiv * zdiv
351               
352               ! tension at T points
353               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
354                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
355                  &   ) * r1_e1e2t(ji,jj)
356               zdt2 = zdt * zdt
357               
358               ! delta at T points
359               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
360
361               ! P/delta at T points
362               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
363               
364               ! stress at T points
365               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
366               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
367             
368            END DO
369         END DO
370         CALL lbc_lnk( zp_delt, 'T', 1. )
371
372         DO jj = 1, jpjm1
373            DO ji = 1, jpim1
374
375               ! P/delta at F points
376               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
377               
378               ! stress at F points
379               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
380
381            END DO
382         END DO
383         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
384 
385
386         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
387         DO jj = 2, jpjm1
388            DO ji = fs_2, fs_jpim1               
389
390               ! U points
391               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
392                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
393                  &                    ) * r1_e2u(ji,jj)                                                                      &
394                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
395                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
396                  &                  ) * r1_e1e2u(ji,jj)
397
398               ! V points
399               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
400                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
401                  &                    ) * r1_e1v(ji,jj)                                                                      &
402                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
403                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
404                  &                  ) * r1_e1e2v(ji,jj)
405
406               ! u_ice at V point
407               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
408                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
409               
410               ! v_ice at U point
411               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
412                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
413
414            END DO
415         END DO
416         !
417         ! --- Computation of ice velocity --- !
418         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
419         !  Bouillon et al. 2009 (eq 34-35) => stable
420         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
421           
422            DO jj = 2, jpjm1
423               DO ji = fs_2, fs_jpim1
424
425                  ! tau_io/(v_oce - v_ice)
426                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
427                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
428
429                  ! tau_bottom/v_ice
430                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
431                  zTauB = - tau_icebfr(ji,jj) / zvel
432
433                  ! Coriolis at V-points (energy conserving formulation)
434                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
435                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
436                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
437
438                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
439                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
440
441                  ! landfast switch => 0 = static friction ; 1 = sliding friction
442                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
443                 
444                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
445                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
446                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
447                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
448                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
449                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
450                     &           ) * zmaskV(ji,jj)
451                  ! Bouillon 2013
452                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
453                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
454                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
455
456               END DO
457            END DO
458            CALL lbc_lnk( v_ice, 'V', -1. )
459           
460#if defined key_agrif
461!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
462            CALL agrif_interp_lim3( 'V' )
463#endif
464            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' )
465
466            DO jj = 2, jpjm1
467               DO ji = fs_2, fs_jpim1
468                               
469                  ! tau_io/(u_oce - u_ice)
470                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
471                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
472
473                  ! tau_bottom/u_ice
474                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
475                  zTauB = - tau_icebfr(ji,jj) / zvel
476
477                  ! Coriolis at U-points (energy conserving formulation)
478                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
479                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
480                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
481                 
482                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
483                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
484
485                  ! landfast switch => 0 = static friction ; 1 = sliding friction
486                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
487
488                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
489                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
490                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
491                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
492                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
493                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
494                     &           ) * zmaskU(ji,jj)
495                  ! Bouillon 2013
496                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
497                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
498                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
499               END DO
500            END DO
501            CALL lbc_lnk( u_ice, 'U', -1. )
502           
503#if defined key_agrif
504!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
505            CALL agrif_interp_lim3( 'U' )
506#endif
507            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' )
508
509         ELSE ! odd iterations
510
511            DO jj = 2, jpjm1
512               DO ji = fs_2, fs_jpim1
513
514                  ! tau_io/(u_oce - u_ice)
515                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
516                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
517
518                  ! tau_bottom/u_ice
519                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
520                  zTauB = - tau_icebfr(ji,jj) / zvel
521
522                  ! Coriolis at U-points (energy conserving formulation)
523                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
524                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
525                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
526
527                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
528                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
529
530                  ! landfast switch => 0 = static friction ; 1 = sliding friction
531                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
532
533                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
534                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
535                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
536                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
537                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
538                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
539                     &           ) * zmaskU(ji,jj)
540                  ! Bouillon 2013
541                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
542                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
543                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
544               END DO
545            END DO
546            CALL lbc_lnk( u_ice, 'U', -1. )
547           
548#if defined key_agrif
549!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
550            CALL agrif_interp_lim3( 'U' )
551#endif
552            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' )
553
554            DO jj = 2, jpjm1
555               DO ji = fs_2, fs_jpim1
556         
557                  ! tau_io/(v_oce - v_ice)
558                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
559                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
560
561                  ! tau_bottom/v_ice
562                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
563                  ztauB = - tau_icebfr(ji,jj) / zvel
564                 
565                  ! Coriolis at V-points (energy conserving formulation)
566                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
567                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
568                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
569
570                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
571                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
572
573                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
574                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
575                 
576                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
577                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
578                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
579                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
580                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
581                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
582                     &           ) * zmaskV(ji,jj)
583                  ! Bouillon 2013
584                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
585                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
586                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
587               END DO
588            END DO
589            CALL lbc_lnk( v_ice, 'V', -1. )
590           
591#if defined key_agrif
592!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
593            CALL agrif_interp_lim3( 'V' )
594#endif
595            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' )
596
597         ENDIF
598         
599         IF(ln_ctl) THEN   ! Convergence test
600            DO jj = 2 , jpjm1
601               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
602            END DO
603            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
604            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
605         ENDIF
606         !
607         !                                                ! ==================== !
608      END DO                                              !  end loop over jter  !
609      !                                                   ! ==================== !
610      !
611      !------------------------------------------------------------------------------!
612      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
613      !------------------------------------------------------------------------------!
614      DO jj = 1, jpjm1
615         DO ji = 1, jpim1
616
617            ! shear at F points
618            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
619               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
620               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
621
622         END DO
623      END DO           
624      CALL lbc_lnk( zds, 'F', 1. )
625     
626      DO jj = 2, jpjm1
627         DO ji = 2, jpim1 ! no vector loop
628           
629            ! tension**2 at T points
630            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
631               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
632               &   ) * r1_e1e2t(ji,jj)
633            zdt2 = zdt * zdt
634           
635            ! shear**2 at T points (doc eq. A16)
636            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
637               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
638               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
639           
640            ! shear at T points
641            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
642
643            ! divergence at T points
644            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
645               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
646               &            ) * r1_e1e2t(ji,jj)
647           
648            ! delta at T points
649            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
650            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
651            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
652
653         END DO
654      END DO
655      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
656     
657      ! --- Store the stress tensor for the next time step --- !
658      stress1_i (:,:) = zs1 (:,:)
659      stress2_i (:,:) = zs2 (:,:)
660      stress12_i(:,:) = zs12(:,:)
661      !
662
663      !------------------------------------------------------------------------------!
664      ! 5) Control prints of residual and charge ellipse
665      !------------------------------------------------------------------------------!
666      !
667      ! print the residual for convergence
668      IF(ln_ctl) THEN
669         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
670         CALL prt_ctl_info(charout)
671         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
672      ENDIF
673
674      ! print charge ellipse
675      ! This can be desactivated once the user is sure that the stress state
676      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
677      IF(ln_ctl) THEN
678         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
679         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
680         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
681         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
682            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
683            CALL prt_ctl_info(charout)
684            DO jj = 2, jpjm1
685               DO ji = 2, jpim1
686                  IF (strength(ji,jj) > 1.0) THEN
687                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
688                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )
689                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
690                     CALL prt_ctl_info(charout)
691                  ENDIF
692               END DO
693            END DO
694            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
695            CALL prt_ctl_info(charout)
696         ENDIF
697      ENDIF
698      !     
699     
700      CALL wrk_dealloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt )
701      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
702      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
703      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
704      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
705
706   END SUBROUTINE lim_rhg
707
708#else
709   !!----------------------------------------------------------------------
710   !!   Default option          Dummy module           NO LIM sea-ice model
711   !!----------------------------------------------------------------------
712CONTAINS
713   SUBROUTINE lim_rhg         ! Dummy routine
714      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?'
715   END SUBROUTINE lim_rhg
716#endif
717
718   !!==============================================================================
719END MODULE limrhg
Note: See TracBrowser for help on using the repository browser.