source: NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE/icedyn_rhg_evp.F90 @ 10888

Last change on this file since 10888 was 10888, checked in by davestorkey, 19 months ago

branches/UKMO/NEMO_4.0_mirror : clear SVN keywords

File size: 58.2 KB
Line 
1MODULE icedyn_rhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_evp  ***
4   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model
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 evp case
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 + mEVP (Bouillon 2013)
13   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
14   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
15   !!----------------------------------------------------------------------
16#if defined key_si3
17   !!----------------------------------------------------------------------
18   !!   'key_si3'                                       SI3 sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
21   !!   rhg_evp_rst     : read/write EVP fields in ice restart
22   !!----------------------------------------------------------------------
23   USE phycst         ! Physical constant
24   USE dom_oce        ! Ocean domain
25   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
26   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
27   USE ice            ! sea-ice: ice variables
28   USE icevar         ! ice_var_sshdyn
29   USE icedyn_rdgrft  ! sea-ice: ice strength
30   USE bdy_oce , ONLY : ln_bdy 
31   USE bdyice 
32#if defined key_agrif
33   USE agrif_ice_interp
34#endif
35   !
36   USE in_out_manager ! I/O manager
37   USE iom            ! I/O manager library
38   USE lib_mpp        ! MPP library
39   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
40   USE lbclnk         ! lateral boundary conditions (or mpp links)
41   USE prtctl         ! Print control
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
47   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
53   !! $Id$
54   !! Software governed by the CeCILL license (see ./LICENSE)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
59      !!-------------------------------------------------------------------
60      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
61      !!                             EVP-C-grid
62      !!
63      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
64      !!  stress and sea-surface slope. Ice-ice interaction is described by
65      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
66      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
67      !!
68      !!  The points in the C-grid look like this, dear reader
69      !!
70      !!                              (ji,jj)
71      !!                                 |
72      !!                                 |
73      !!                      (ji-1,jj)  |  (ji,jj)
74      !!                             ---------   
75      !!                            |         |
76      !!                            | (ji,jj) |------(ji,jj)
77      !!                            |         |
78      !!                             ---------   
79      !!                     (ji-1,jj-1)     (ji,jj-1)
80      !!
81      !! ** Inputs  : - wind forcing (stress), oceanic currents
82      !!                ice total volume (vt_i) per unit area
83      !!                snow total volume (vt_s) per unit area
84      !!
85      !! ** Action  : - compute u_ice, v_ice : the components of the
86      !!                sea-ice velocity vector
87      !!              - compute delta_i, shear_i, divu_i, which are inputs
88      !!                of the ice thickness distribution
89      !!
90      !! ** Steps   : 0) compute mask at F point
91      !!              1) Compute ice snow mass, ice strength
92      !!              2) Compute wind, oceanic stresses, mass terms and
93      !!                 coriolis terms of the momentum equation
94      !!              3) Solve the momentum equation (iterative procedure)
95      !!              4) Recompute delta, shear and divergence
96      !!                 (which are inputs of the ITD) & store stress
97      !!                 for the next time step
98      !!              5) Diagnostics including charge ellipse
99      !!
100      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
101      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
102      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
103      !!              (i.e. more stable and better convergence)
104      !!
105      !! References : Hunke and Dukowicz, JPO97
106      !!              Bouillon et al., Ocean Modelling 2009
107      !!              Bouillon et al., Ocean Modelling 2013
108      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
109      !!-------------------------------------------------------------------
110      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
111      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
113      !!
114      LOGICAL, PARAMETER ::   ll_bdy_substep = .FALSE. ! temporary option to call bdy at each sub-time step (T)
115      !                                                                              or only at the main time step (F)
116      INTEGER ::   ji, jj       ! dummy loop indices
117      INTEGER ::   jter         ! local integers
118      !
119      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
120      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
121      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
122      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
123      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
124      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
125      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars
126      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
127      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
128      !
129      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
130      REAL(wp) ::   zintb, zintn                                        ! dummy argument
131      REAL(wp) ::   zfac_x, zfac_y
132      REAL(wp) ::   zshear, zdum1, zdum2
133      !
134      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
135      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
136      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
137      !
138      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
139      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
140      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
141      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
142      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
143      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param)
144      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
145      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
146      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
147      !
148      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
149      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
150!!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
151      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
152      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
153      !                                                                 !    ice bottom surface if ice is embedded   
154      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
155      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
156      !
157      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
158      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
159      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
160
161      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
162      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
163      !! --- diags
164      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
166      !! --- SIMIP diags
167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress 
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
179      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
180      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
181      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
182      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
183      !!-------------------------------------------------------------------
184
185      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
186      !
187!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
188      !------------------------------------------------------------------------------!
189      ! 0) mask at F points for the ice
190      !------------------------------------------------------------------------------!
191      ! ocean/land mask
192      DO jj = 1, jpjm1
193         DO ji = 1, jpim1      ! NO vector opt.
194            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
195         END DO
196      END DO
197      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
198
199      ! Lateral boundary conditions on velocity (modify zfmask)
200      zwf(:,:) = zfmask(:,:)
201      DO jj = 2, jpjm1
202         DO ji = fs_2, fs_jpim1   ! vector opt.
203            IF( zfmask(ji,jj) == 0._wp ) THEN
204               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
205            ENDIF
206         END DO
207      END DO
208      DO jj = 2, jpjm1
209         IF( zfmask(1,jj) == 0._wp ) THEN
210            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
211         ENDIF
212         IF( zfmask(jpi,jj) == 0._wp ) THEN
213            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
214         ENDIF
215      END DO
216      DO ji = 2, jpim1
217         IF( zfmask(ji,1) == 0._wp ) THEN
218            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
219         ENDIF
220         IF( zfmask(ji,jpj) == 0._wp ) THEN
221            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
222         ENDIF
223      END DO
224      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
225
226      !------------------------------------------------------------------------------!
227      ! 1) define some variables and initialize arrays
228      !------------------------------------------------------------------------------!
229      zrhoco = rau0 * rn_cio 
230
231      ! ecc2: square of yield ellipse eccenticrity
232      ecc2    = rn_ecc * rn_ecc
233      z1_ecc2 = 1._wp / ecc2
234
235      ! Time step for subcycling
236      zdtevp   = rdt_ice / REAL( nn_nevp )
237      z1_dtevp = 1._wp / zdtevp
238
239      ! alpha parameters (Bouillon 2009)
240      IF( .NOT. ln_aEVP ) THEN
241         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
242         zalph2 = zalph1 * z1_ecc2
243
244         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
245         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
246      ENDIF
247         
248      ! Initialise stress tensor
249      zs1 (:,:) = pstress1_i (:,:) 
250      zs2 (:,:) = pstress2_i (:,:)
251      zs12(:,:) = pstress12_i(:,:)
252
253      ! Ice strength
254      CALL ice_strength
255
256      ! scale factors
257      DO jj = 2, jpjm1
258         DO ji = fs_2, fs_jpim1
259            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
260            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
261         END DO
262      END DO
263
264      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
265      IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile
266      ELSE                                               ;   zkt = 0._wp
267      ENDIF
268      !
269      !------------------------------------------------------------------------------!
270      ! 2) Wind / ocean stress, mass terms, coriolis terms
271      !------------------------------------------------------------------------------!
272      ! sea surface height
273      !    embedded sea ice: compute representative ice top surface
274      !    non-embedded sea ice: use ocean surface for slope calculation
275      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
276
277      DO jj = 2, jpjm1
278         DO ji = fs_2, fs_jpim1
279
280            ! ice fraction at U-V points
281            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)
282            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)
283
284            ! Ice/snow mass at U-V points
285            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
286            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
287            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
288            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
289            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
290
291            ! Ocean currents at U-V points
292            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
293               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
294           
295            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
296               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
297
298            ! Coriolis at T points (m*f)
299            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
300
301            ! dt/m at T points (for alpha and beta coefficients)
302            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
303           
304            ! m/dt
305            zmU_t(ji,jj)    = zmassU * z1_dtevp
306            zmV_t(ji,jj)    = zmassV * z1_dtevp
307           
308            ! Drag ice-atm.
309            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
310            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
311
312            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
313            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
314            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
315
316            ! masks
317            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
318            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
319
320            ! switches
321            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
322            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
323
324         END DO
325      END DO
326      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. )
327      !
328      !                                  !== Landfast ice parameterization ==!
329      !
330      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
331         DO jj = 2, jpjm1
332            DO ji = fs_2, fs_jpim1
333               ! ice thickness at U-V points
334               zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
335               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
336               ! ice-bottom stress at U points
337               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
338               zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
339               ! ice-bottom stress at V points
340               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
341               zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
342               ! ice_bottom stress at T points
343               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
344               tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
345            END DO
346         END DO
347         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. )
348         !
349      ELSEIF( ln_landfast_home ) THEN          !-- Home made
350         DO jj = 2, jpjm1
351            DO ji = fs_2, fs_jpim1
352               zTauU_ib(ji,jj) = tau_icebfr(ji,jj)
353               zTauV_ib(ji,jj) = tau_icebfr(ji,jj)
354            END DO
355         END DO
356         !
357      ELSE                                     !-- no landfast
358         DO jj = 2, jpjm1
359            DO ji = fs_2, fs_jpim1
360               zTauU_ib(ji,jj) = 0._wp
361               zTauV_ib(ji,jj) = 0._wp
362            END DO
363         END DO
364      ENDIF
365      IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )
366
367      !------------------------------------------------------------------------------!
368      ! 3) Solution of the momentum equation, iterative procedure
369      !------------------------------------------------------------------------------!
370      !
371      !                                               !----------------------!
372      DO jter = 1 , nn_nevp                           !    loop over jter    !
373         !                                            !----------------------!       
374         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
375         !
376!!$         IF(ln_ctl) THEN   ! Convergence test
377!!$            DO jj = 1, jpjm1
378!!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
379!!$               zv_ice(:,jj) = v_ice(:,jj)
380!!$            END DO
381!!$         ENDIF
382
383         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
384         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
385            DO ji = 1, jpim1
386
387               ! shear at F points
388               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)   &
389                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
390                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
391
392            END DO
393         END DO
394         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
395
396         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
397            DO ji = 2, jpi ! no vector loop
398
399               ! shear**2 at T points (doc eq. A16)
400               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
401                  &   + 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)  &
402                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
403             
404               ! divergence at T points
405               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
406                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
407                  &    ) * r1_e1e2t(ji,jj)
408               zdiv2 = zdiv * zdiv
409               
410               ! tension at T points
411               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)   &
412                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
413                  &   ) * r1_e1e2t(ji,jj)
414               zdt2 = zdt * zdt
415               
416               ! delta at T points
417               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
418
419               ! P/delta at T points
420               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
421
422               ! alpha & beta for aEVP
423               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
424               !   alpha = beta = sqrt(4*gamma)
425               IF( ln_aEVP ) THEN
426                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
427                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
428                  zalph2   = zalph1
429                  z1_alph2 = z1_alph1
430               ENDIF
431               
432               ! stress at T points (zkt/=0 if landfast)
433               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
434               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
435             
436            END DO
437         END DO
438         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
439
440         DO jj = 1, jpjm1
441            DO ji = 1, jpim1
442
443               ! alpha & beta for aEVP
444               IF( ln_aEVP ) THEN
445                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
446                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
447                  zbeta(ji,jj) = zalph2
448               ENDIF
449               
450               ! P/delta at F points
451               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) )
452               
453               ! stress at F points (zkt/=0 if landfast)
454               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
455
456            END DO
457         END DO
458
459         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
460         DO jj = 2, jpjm1
461            DO ji = fs_2, fs_jpim1               
462               !                   !--- U points
463               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
464                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
465                  &                    ) * r1_e2u(ji,jj)                                                                      &
466                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
467                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
468                  &                  ) * r1_e1e2u(ji,jj)
469               !
470               !                !--- V points
471               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
472                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
473                  &                    ) * r1_e1v(ji,jj)                                                                      &
474                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
475                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
476                  &                  ) * r1_e1e2v(ji,jj)
477               !
478               !                !--- u_ice at V point
479               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
480                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
481               !
482               !                !--- v_ice at U point
483               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
484                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
485            END DO
486         END DO
487         !
488         ! --- Computation of ice velocity --- !
489         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
490         !  Bouillon et al. 2009 (eq 34-35) => stable
491         IF( MOD(jter,2) == 0 ) THEN ! even iterations
492            !
493            DO jj = 2, jpjm1
494               DO ji = fs_2, fs_jpim1
495                  !                 !--- tau_io/(v_oce - v_ice)
496                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
497                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
498                  !                 !--- Ocean-to-Ice stress
499                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
500                  !
501                  !                 !--- tau_bottom/v_ice
502                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
503                  zTauB = - zTauV_ib(ji,jj) / zvel
504                  !
505                  !                 !--- Coriolis at V-points (energy conserving formulation)
506                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
507                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
508                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
509                  !
510                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
511                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
512                  !
513                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
514                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
515                  !
516                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
517                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
518                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
519                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
520                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
521                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
522                     &           ) * zmaskV(ji,jj)
523                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
524                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
525                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
526                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
527                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
528                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
529                     &            ) * zmaskV(ji,jj)
530                  ENDIF
531               END DO
532            END DO
533            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
534            !
535#if defined key_agrif
536!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
537            CALL agrif_interp_ice( 'V' )
538#endif
539            IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' )
540            !
541            DO jj = 2, jpjm1
542               DO ji = fs_2, fs_jpim1         
543                  !                 !--- tau_io/(u_oce - u_ice)
544                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
545                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
546                  !                 !--- Ocean-to-Ice stress
547                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
548                  !
549                  !                 !--- tau_bottom/u_ice
550                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
551                  zTauB = - zTauU_ib(ji,jj) / zvel
552                  !
553                  !                 !--- Coriolis at U-points (energy conserving formulation)
554                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
555                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
556                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
557                  !
558                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
559                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
560                  !
561                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
562                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
563                  !
564                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
565                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
566                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
567                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
568                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
569                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
570                     &            ) * zmaskU(ji,jj)
571                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
572                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
573                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
574                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
575                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
576                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
577                     &            ) * zmaskU(ji,jj)
578                  ENDIF
579               END DO
580            END DO
581            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
582            !
583#if defined key_agrif
584!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
585            CALL agrif_interp_ice( 'U' )
586#endif
587            IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' )
588            !
589         ELSE ! odd iterations
590            !
591            DO jj = 2, jpjm1
592               DO ji = fs_2, fs_jpim1
593                  !                 !--- tau_io/(u_oce - u_ice)
594                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
595                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
596                  !                 !--- Ocean-to-Ice stress
597                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
598                  !
599                  !                 !--- tau_bottom/u_ice
600                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
601                  zTauB = - zTauU_ib(ji,jj) / zvel
602                  !
603                  !                 !--- Coriolis at U-points (energy conserving formulation)
604                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
605                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
606                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
607                  !
608                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
609                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
610                  !
611                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
612                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
613                  !
614                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
615                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
616                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
617                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
618                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
619                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
620                     &            ) * zmaskU(ji,jj)
621                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
622                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
623                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
624                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
625                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
626                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
627                     &            ) * zmaskU(ji,jj)
628                  ENDIF
629               END DO
630            END DO
631            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
632            !
633#if defined key_agrif
634!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
635            CALL agrif_interp_ice( 'U' )
636#endif
637            IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' )
638            !
639            DO jj = 2, jpjm1
640               DO ji = fs_2, fs_jpim1
641                  !                 !--- tau_io/(v_oce - v_ice)
642                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
643                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
644                  !                 !--- Ocean-to-Ice stress
645                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
646                  !
647                  !                 !--- tau_bottom/v_ice
648                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
649                  zTauB = - zTauV_ib(ji,jj) / zvel
650                  !
651                  !                 !--- Coriolis at v-points (energy conserving formulation)
652                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
653                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
654                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
655                  !
656                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
657                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
658                  !
659                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
660                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
661                  !
662                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
663                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
664                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
665                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
666                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
667                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
668                     &           ) * zmaskV(ji,jj)
669                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
670                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
671                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
672                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
673                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
674                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
675                     &            ) * zmaskV(ji,jj)
676                  ENDIF
677               END DO
678            END DO
679            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
680            !
681#if defined key_agrif
682!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
683            CALL agrif_interp_ice( 'V' )
684#endif
685            IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' )
686            !
687         ENDIF
688
689!!$         IF(ln_ctl) THEN   ! Convergence test
690!!$            DO jj = 2 , jpjm1
691!!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
692!!$            END DO
693!!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
694!!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
695!!$         ENDIF
696         !
697         !                                                ! ==================== !
698      END DO                                              !  end loop over jter  !
699      !                                                   ! ==================== !
700      !
701      IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN
702         CALL bdy_ice_dyn( 'U' )
703         CALL bdy_ice_dyn( 'V' )
704      ENDIF
705      !
706      !------------------------------------------------------------------------------!
707      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
708      !------------------------------------------------------------------------------!
709      DO jj = 1, jpjm1
710         DO ji = 1, jpim1
711
712            ! shear at F points
713            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)   &
714               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
715               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
716
717         END DO
718      END DO           
719     
720      DO jj = 2, jpjm1
721         DO ji = 2, jpim1 ! no vector loop
722           
723            ! tension**2 at T points
724            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)   &
725               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
726               &   ) * r1_e1e2t(ji,jj)
727            zdt2 = zdt * zdt
728           
729            ! shear**2 at T points (doc eq. A16)
730            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
731               &   + 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)  &
732               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
733           
734            ! shear at T points
735            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
736
737            ! divergence at T points
738            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
739               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
740               &             ) * r1_e1e2t(ji,jj)
741           
742            ! delta at T points
743            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
744            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
745            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
746
747         END DO
748      END DO
749      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
750     
751      ! --- Store the stress tensor for the next time step --- !
752      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
753      pstress1_i (:,:) = zs1 (:,:)
754      pstress2_i (:,:) = zs2 (:,:)
755      pstress12_i(:,:) = zs12(:,:)
756      !
757
758      !------------------------------------------------------------------------------!
759      ! 5) diagnostics
760      !------------------------------------------------------------------------------!
761      DO jj = 1, jpj
762         DO ji = 1, jpi
763            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
764         END DO
765      END DO
766
767      ! --- divergence, shear and strength --- !
768      IF( iom_use('icediv') )   CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) )   ! divergence
769      IF( iom_use('iceshe') )   CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) )   ! shear
770      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
771
772      ! --- charge ellipse --- !
773      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
774         !
775         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
776         !         
777         DO jj = 2, jpjm1
778            DO ji = 2, jpim1
779               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
780                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
781                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
782
783               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
784
785               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
786
787!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
788!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
789!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
790!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
791               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
792               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
793               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
794            END DO
795         END DO
796         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
797         !
798         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
799         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
800         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
801         !
802         DEALLOCATE( zsig1 , zsig2 , zsig3 )
803      ENDIF
804     
805      ! --- SIMIP --- !
806      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
807         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
808         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
809         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
810
811         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
812            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
813            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
814            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
815         
816         DO jj = 2, jpjm1
817            DO ji = 2, jpim1
818               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
819               
820               ! Stress tensor invariants (normal and shear stress N/m)
821               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
822               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
823               
824               ! Stress terms of the momentum equation (N/m2)
825               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term
826               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
827               
828               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term
829               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
830               
831               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term
832               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
833               
834               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
835               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
836               
837               ! 2D ice mass, snow mass, area transport arrays (X, Y)
838               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
839               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
840               
841               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
842               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
843               
844               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
845               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
846               
847               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
848               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
849               
850            END DO
851         END DO
852         
853         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
854            &                zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
855            &                zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
856            &                zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
857                 
858         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi  , 'U', -1., zdiag_vtau_oi  , 'V', -1.,   &
859            &                zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
860            &                zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
861            &                zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
862         
863         IF( iom_use('normstr' ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
864         IF( iom_use('sheastr' ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
865         IF( iom_use('dssh_dx' ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
866         IF( iom_use('dssh_dy' ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
867         IF( iom_use('corstrx' ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
868         IF( iom_use('corstry' ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
869         IF( iom_use('intstrx' ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
870         IF( iom_use('intstry' ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
871         IF( iom_use('utau_oi' ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
872         IF( iom_use('vtau_oi' ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
873         IF( iom_use('xmtrpice') )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
874         IF( iom_use('ymtrpice') )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport
875         IF( iom_use('xmtrpsnw') )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
876         IF( iom_use('ymtrpsnw') )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
877         IF( iom_use('xatrp'   ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
878         IF( iom_use('yatrp'   ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
879
880         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
881            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
882            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
883            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
884
885      ENDIF
886      !
887   END SUBROUTINE ice_dyn_rhg_evp
888
889
890   SUBROUTINE rhg_evp_rst( cdrw, kt )
891      !!---------------------------------------------------------------------
892      !!                   ***  ROUTINE rhg_evp_rst  ***
893      !!                     
894      !! ** Purpose :   Read or write RHG file in restart file
895      !!
896      !! ** Method  :   use of IOM library
897      !!----------------------------------------------------------------------
898      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
899      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
900      !
901      INTEGER  ::   iter            ! local integer
902      INTEGER  ::   id1, id2, id3   ! local integers
903      !!----------------------------------------------------------------------
904      !
905      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
906         !                                   ! ---------------
907         IF( ln_rstart ) THEN                   !* Read the restart file
908            !
909            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
910            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
911            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
912            !
913            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
914               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
915               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
916               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
917            ELSE                                     ! start rheology from rest
918               IF(lwp) WRITE(numout,*)
919               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
920               stress1_i (:,:) = 0._wp
921               stress2_i (:,:) = 0._wp
922               stress12_i(:,:) = 0._wp
923            ENDIF
924         ELSE                                   !* Start from rest
925            IF(lwp) WRITE(numout,*)
926            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
927            stress1_i (:,:) = 0._wp
928            stress2_i (:,:) = 0._wp
929            stress12_i(:,:) = 0._wp
930         ENDIF
931         !
932      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
933         !                                   ! -------------------
934         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
935         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
936         !
937         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
938         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
939         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
940         !
941      ENDIF
942      !
943   END SUBROUTINE rhg_evp_rst
944
945#else
946   !!----------------------------------------------------------------------
947   !!   Default option         Empty module           NO SI3 sea-ice model
948   !!----------------------------------------------------------------------
949#endif
950
951   !!==============================================================================
952END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.