New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limrhg.F90 in branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 7293

Last change on this file since 7293 was 7158, checked in by clem, 7 years ago

debug branch

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