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.
icedyn_rhg_evp.F90 in NEMO/branches/UKMO/NEMO_4.0_mirror/src/ICE – NEMO

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

Last change on this file since 11081 was 11081, checked in by davestorkey, 5 years ago

UKMO/NEMO_4.0_mirror : update to be a copy of rev 11079 of release-4.0.

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