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.1_fix_cpl_v2/src/ICE – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icedyn_rhg_evp.F90 @ 12937

Last change on this file since 12937 was 12937, checked in by dancopsey, 3 years ago

Merge in Clem's branch. It was originally here:

svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/UKMO/NEMO_4.0.1_dan_test_clems_branch

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