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/2020/SI3-03_VP_rheology/src/ICE – NEMO

source: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90 @ 13536

Last change on this file since 13536 was 13536, checked in by vancop, 20 months ago

VP rheology, first code version finalised

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