source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg_evp.F90 @ 13279

Last change on this file since 13279 was 13279, checked in by clem, 3 months ago

merge with r4.0-HEAD at r13278

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