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/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_rhg_evp.F90 @ 15493

Last change on this file since 15493 was 15493, checked in by clem, 8 months ago

4.0-HEAD: correct output for rheology

  • Property svn:keywords set to Id
File size: 63.9 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) ::   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) ::   zdelta, zp_delt                 ! delta and P/delta at T points
139      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension
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, zfac, z1_strength
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p         
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( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) 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
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
387         DO jj = 2, jpjm1
388            DO ji = 2, jpim1 ! no vector loop
389
390               ! shear**2 at T points (doc eq. A16)
391               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
392                  &   + 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)  &
393                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
394             
395               ! divergence at T points
396               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
397                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
398                  &    ) * r1_e1e2t(ji,jj)
399               zdiv2 = zdiv * zdiv
400               
401               ! tension at T points
402               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)   &
403                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
404                  &   ) * r1_e1e2t(ji,jj)
405               zdt2 = zdt * zdt
406               
407               ! delta at T points
408               zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
409
410            END DO
411         END DO
412         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1._wp )
413         
414         ! P/delta at T points
415         DO jj = 1, jpj
416            DO ji = 1, jpi
417               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
418            END DO
419         END DO
420
421         DO jj = 2, jpj    ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
422            DO ji = 2, jpi ! no vector loop
423
424               ! divergence at T points (duplication to avoid communications)
425               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
426                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
427                  &    ) * r1_e1e2t(ji,jj)
428               
429               ! tension at T points (duplication to avoid communications)
430               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)   &
431                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
432                  &   ) * r1_e1e2t(ji,jj)
433
434               ! alpha for aEVP
435               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
436               !   alpha = beta = sqrt(4*gamma)
437               IF( ln_aEVP ) THEN
438                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
439                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
440                  zalph2   = zalph1
441                  z1_alph2 = z1_alph1
442                  ! explicit:
443                  ! z1_alph1 = 1._wp / zalph1
444                  ! z1_alph2 = 1._wp / zalph1
445                  ! zalph1 = zalph1 - 1._wp
446                  ! zalph2 = zalph1
447               ENDIF
448               
449               ! stress at T points (zkt/=0 if landfast)
450               zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1
451               zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
452             
453            END DO
454         END DO
455
456         ! Save beta at T-points for further computations
457         IF( ln_aEVP ) THEN
458            DO jj = 1, jpj
459               DO ji = 1, jpi
460                  zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
461               END DO
462            END DO
463         ENDIF
464         
465         DO jj = 1, jpjm1
466            DO ji = 1, jpim1
467
468               ! alpha for aEVP
469               IF( ln_aEVP ) THEN
470                  zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) )
471                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
472                  ! explicit:
473                  ! z1_alph2 = 1._wp / zalph2
474                  ! zalph2 = zalph2 - 1._wp
475               ENDIF
476               
477               ! P/delta at F points
478               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) )
479               
480               ! stress at F points (zkt/=0 if landfast)
481               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
482
483            END DO
484         END DO
485
486         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
487         DO jj = 2, jpjm1
488            DO ji = fs_2, fs_jpim1               
489               !                   !--- U points
490               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
491                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
492                  &                    ) * r1_e2u(ji,jj)                                                                      &
493                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
494                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
495                  &                  ) * r1_e1e2u(ji,jj)
496               !
497               !                !--- V points
498               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
499                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
500                  &                    ) * r1_e1v(ji,jj)                                                                      &
501                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
502                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
503                  &                  ) * r1_e1e2v(ji,jj)
504               !
505               !                !--- ice currents at U-V point
506               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)
507               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)
508               !
509            END DO
510         END DO
511         !
512         ! --- Computation of ice velocity --- !
513         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
514         !  Bouillon et al. 2009 (eq 34-35) => stable
515         IF( MOD(jter,2) == 0 ) THEN ! even iterations
516            !
517            DO jj = 2, jpjm1
518               DO ji = fs_2, fs_jpim1
519                  !                 !--- tau_io/(v_oce - v_ice)
520                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
521                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
522                  !                 !--- Ocean-to-Ice stress
523                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
524                  !
525                  !                 !--- tau_bottom/v_ice
526                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
527                  zTauB = ztauy_base(ji,jj) / zvel
528                  !                 !--- OceanBottom-to-Ice stress
529                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
530                  !
531                  !                 !--- Coriolis at V-points (energy conserving formulation)
532                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
533                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
534                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
535                  !
536                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
537                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
538                  !
539                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
540                  !                                         1 = sliding friction : TauB < RHS
541                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
542                  !
543                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
544                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
545                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
546                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
547                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
548                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
549                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
550                        &                                    ) / ( zbetav + 1._wp )                                              &
551                        &             ) * 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
552                        &           )   * zmsk00y(ji,jj)
553                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
554                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
555                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
556                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
557                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
558                        &             ) * 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
559                        &            )  * zmsk00y(ji,jj)
560                  ENDIF
561               END DO
562            END DO
563            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
564            !
565#if defined key_agrif
566!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
567            CALL agrif_interp_ice( 'V' )
568#endif
569            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
570            !
571            DO jj = 2, jpjm1
572               DO ji = fs_2, fs_jpim1         
573                  !                 !--- tau_io/(u_oce - u_ice)
574                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
575                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
576                  !                 !--- Ocean-to-Ice stress
577                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
578                  !
579                  !                 !--- tau_bottom/u_ice
580                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
581                  zTauB = ztaux_base(ji,jj) / zvel
582                  !                 !--- OceanBottom-to-Ice stress
583                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
584                  !
585                  !                 !--- Coriolis at U-points (energy conserving formulation)
586                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
587                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
588                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
589                  !
590                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
591                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
592                  !
593                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
594                  !                                         1 = sliding friction : TauB < RHS
595                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
596                  !
597                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
598                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
599                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
600                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
601                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
602                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
603                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
604                        &                                    ) / ( zbetau + 1._wp )                                              &
605                        &             ) * 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
606                        &           )   * zmsk00x(ji,jj)
607                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
608                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
609                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
610                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
611                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
612                        &             ) * 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
613                        &           )   * zmsk00x(ji,jj)
614                  ENDIF
615               END DO
616            END DO
617            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
618            !
619#if defined key_agrif
620!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
621            CALL agrif_interp_ice( 'U' )
622#endif
623            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
624            !
625         ELSE ! odd iterations
626            !
627            DO jj = 2, jpjm1
628               DO ji = fs_2, fs_jpim1
629                  !                 !--- tau_io/(u_oce - u_ice)
630                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
631                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
632                  !                 !--- Ocean-to-Ice stress
633                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
634                  !
635                  !                 !--- tau_bottom/u_ice
636                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
637                  zTauB = ztaux_base(ji,jj) / zvel
638                  !                 !--- OceanBottom-to-Ice stress
639                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
640                  !
641                  !                 !--- Coriolis at U-points (energy conserving formulation)
642                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
643                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
644                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
645                  !
646                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
647                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
648                  !
649                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
650                  !                                         1 = sliding friction : TauB < RHS
651                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
652                  !
653                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
654                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
655                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
656                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
657                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
658                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
659                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
660                        &                                    ) / ( zbetau + 1._wp )                                              &
661                        &             ) * 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
662                        &           )   * zmsk00x(ji,jj)
663                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
664                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
665                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
666                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
667                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
668                        &             ) * 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
669                        &           )   * zmsk00x(ji,jj)
670                  ENDIF
671               END DO
672            END DO
673            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
674            !
675#if defined key_agrif
676!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
677            CALL agrif_interp_ice( 'U' )
678#endif
679            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
680            !
681            DO jj = 2, jpjm1
682               DO ji = fs_2, fs_jpim1
683                  !                 !--- tau_io/(v_oce - v_ice)
684                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
685                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
686                  !                 !--- Ocean-to-Ice stress
687                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
688                  !
689                  !                 !--- tau_bottom/v_ice
690                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
691                  zTauB = ztauy_base(ji,jj) / zvel
692                  !                 !--- OceanBottom-to-Ice stress
693                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
694                  !
695                  !                 !--- Coriolis at v-points (energy conserving formulation)
696                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
697                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
698                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
699                  !
700                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
701                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
702                  !
703                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
704                  !                                         1 = sliding friction : TauB < RHS
705                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
706                  !
707                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
708                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
709                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(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) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
712                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
713                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
714                        &                                    ) / ( zbetav + 1._wp )                                              & 
715                        &             ) * 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
716                        &           )   * zmsk00y(ji,jj)
717                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
718                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
719                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
720                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
721                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
722                        &             ) * 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
723                        &           )   * zmsk00y(ji,jj)
724                  ENDIF
725               END DO
726            END DO
727            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
728            !
729#if defined key_agrif
730!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
731            CALL agrif_interp_ice( 'V' )
732#endif
733            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
734            !
735         ENDIF
736
737         ! convergence test
738         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
739         !
740         !                                                ! ==================== !
741      END DO                                              !  end loop over jter  !
742      !                                                   ! ==================== !
743      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta )
744      !
745      !------------------------------------------------------------------------------!
746      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
747      !------------------------------------------------------------------------------!
748      DO jj = 1, jpjm1
749         DO ji = 1, jpim1
750
751            ! shear at F points
752            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)   &
753               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
754               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
755
756         END DO
757      END DO           
758     
759      DO jj = 2, jpjm1
760         DO ji = 2, jpim1 ! no vector loop
761           
762            ! tension**2 at T points
763            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)   &
764               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
765               &   ) * r1_e1e2t(ji,jj)
766            zdt2 = zdt * zdt
767           
768            zten_i(ji,jj) = zdt
769
770            ! shear**2 at T points (doc eq. A16)
771            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
772               &   + 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)  &
773               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
774           
775            ! shear at T points
776            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
777
778            ! divergence at T points
779            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
780               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
781               &             ) * r1_e1e2t(ji,jj)
782           
783            ! delta at T points
784            zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 
785            rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
786            pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
787
788         END DO
789      END DO
790      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., &
791         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1. )
792     
793      ! --- Store the stress tensor for the next time step --- !
794      pstress1_i (:,:) = zs1 (:,:)
795      pstress2_i (:,:) = zs2 (:,:)
796      pstress12_i(:,:) = zs12(:,:)
797      !
798
799      !------------------------------------------------------------------------------!
800      ! 5) diagnostics
801      !------------------------------------------------------------------------------!
802      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
803      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
804         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
805         !
806         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
807            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
808         !
809         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
810         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
811         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
812         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
813         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
814         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
815      ENDIF
816       
817      ! --- divergence, shear and strength --- !
818      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
819      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
820      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
821
822      ! --- Stress tensor invariants (SIMIP diags) --- !
823      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
824         !
825         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
826         !         
827         DO jj = 1, jpj
828            DO ji = 1, jpi
829           
830               ! Ice stresses
831               ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
832               ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
833               ! I know, this can be confusing...
834               zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 
835               zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
836               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
837               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
838               
839               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
840               zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                      ! 1st stress invariant, aka average normal stress, aka negative pressure
841               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )  ! 2nd  ''       ''    , aka maximum shear stress
842               
843            END DO
844         END DO         
845         !
846         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
847         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
848         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
849         
850         DEALLOCATE ( zsig_I, zsig_II )
851         
852      ENDIF
853
854      ! --- Normalized stress tensor principal components --- !
855      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
856      ! Recommendation 1 : we use ice strength, not replacement pressure
857      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
858!!$      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
859!!$         !
860!!$         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )         
861!!$         !         
862!!$         DO jj = 1, jpj
863!!$            DO ji = 1, jpi
864!!$           
865!!$               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
866!!$               !                        and **deformations** at current iterates
867!!$               !                        following Lemieux & Dupont (2020)
868!!$               zfac             =   zp_delt(ji,jj)
869!!$               zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
870!!$               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
871!!$               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
872!!$               
873!!$               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
874!!$               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                      ! 1st stress invariant, aka average normal stress, aka negative pressure
875!!$               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )  ! 2nd  ''       ''    , aka maximum shear stress
876!!$     
877!!$               ! Normalized  principal stresses (used to display the ellipse)
878!!$               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
879!!$               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
880!!$               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
881!!$            END DO
882!!$         END DO               
883!!$         !
884!!$         CALL iom_put( 'sig1_pnorm' , zsig1_p )
885!!$         CALL iom_put( 'sig2_pnorm' , zsig2_p )
886!!$     
887!!$         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
888!!$         
889!!$      ENDIF
890
891      ! --- SIMIP --- !
892      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
893         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
894         !
895         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., &
896            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
897
898         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
899         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
900         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
901         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
902         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
903         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
904      ENDIF
905
906      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
907         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
908         !
909         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
910            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
911         !
912         DO jj = 2, jpjm1
913            DO ji = 2, jpim1
914               ! 2D ice mass, snow mass, area transport arrays (X, Y)
915               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
916               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
917
918               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
919               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
920
921               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
922               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
923
924               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
925               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
926
927            END DO
928         END DO
929
930         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
931            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
932            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
933
934         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
935         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport
936         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
937         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
938         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
939         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
940
941         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
942            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
943
944      ENDIF
945      !
946      ! --- convergence tests --- !
947      IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
948         IF( iom_use('uice_cvg') ) THEN
949            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
950               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
951                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
952            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
953               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
954                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
955            ENDIF
956         ENDIF
957      ENDIF     
958      !
959      DEALLOCATE( zmsk00, zmsk15 )
960      !
961   END SUBROUTINE ice_dyn_rhg_evp
962
963
964   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb )
965      !!----------------------------------------------------------------------
966      !!                    ***  ROUTINE rhg_cvg  ***
967      !!                     
968      !! ** Purpose :   check convergence of oce rheology
969      !!
970      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity
971      !!                during the sub timestepping of rheology so as:
972      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) )
973      !!                This routine is called every sub-iteration, so it is cpu expensive
974      !!
975      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)   
976      !!----------------------------------------------------------------------
977      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
978      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
979      !!
980      INTEGER           ::   it, idtime, istatus
981      INTEGER           ::   ji, jj          ! dummy loop indices
982      REAL(wp)          ::   zresm           ! local real
983      CHARACTER(len=20) ::   clname
984      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence
985      !!----------------------------------------------------------------------
986
987      ! create file
988      IF( kt == nit000 .AND. kiter == 1 ) THEN
989         !
990         IF( lwp ) THEN
991            WRITE(numout,*)
992            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
993            WRITE(numout,*) '~~~~~~~'
994         ENDIF
995         !
996         IF( lwm ) THEN
997            clname = 'ice_cvg.nc'
998            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
999            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
1000            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
1001            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid )
1002            istatus = NF90_ENDDEF(ncvgid)
1003         ENDIF
1004         !
1005      ENDIF
1006
1007      ! time
1008      it = ( kt - 1 ) * kitermax + kiter
1009     
1010      ! convergence
1011      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
1012         zresm = 0._wp
1013      ELSE
1014         DO jj = 1, jpj
1015            DO ji = 1, jpi
1016               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
1017                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)
1018            END DO
1019         END DO
1020         zresm = MAXVAL( zres )
1021         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
1022      ENDIF
1023
1024      IF( lwm ) THEN
1025         ! write variables
1026         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
1027         ! close file
1028         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid)
1029      ENDIF
1030     
1031   END SUBROUTINE rhg_cvg
1032
1033
1034   SUBROUTINE rhg_evp_rst( cdrw, kt )
1035      !!---------------------------------------------------------------------
1036      !!                   ***  ROUTINE rhg_evp_rst  ***
1037      !!                     
1038      !! ** Purpose :   Read or write RHG file in restart file
1039      !!
1040      !! ** Method  :   use of IOM library
1041      !!----------------------------------------------------------------------
1042      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1043      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
1044      !
1045      INTEGER  ::   iter            ! local integer
1046      INTEGER  ::   id1, id2, id3   ! local integers
1047      !!----------------------------------------------------------------------
1048      !
1049      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
1050         !                                   ! ---------------
1051         IF( ln_rstart ) THEN                   !* Read the restart file
1052            !
1053            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
1054            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
1055            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
1056            !
1057            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
1058               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
1059               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
1060               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
1061            ELSE                                     ! start rheology from rest
1062               IF(lwp) WRITE(numout,*)
1063               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
1064               stress1_i (:,:) = 0._wp
1065               stress2_i (:,:) = 0._wp
1066               stress12_i(:,:) = 0._wp
1067            ENDIF
1068         ELSE                                   !* Start from rest
1069            IF(lwp) WRITE(numout,*)
1070            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
1071            stress1_i (:,:) = 0._wp
1072            stress2_i (:,:) = 0._wp
1073            stress12_i(:,:) = 0._wp
1074         ENDIF
1075         !
1076      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1077         !                                   ! -------------------
1078         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
1079         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
1080         !
1081         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
1082         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
1083         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
1084         !
1085      ENDIF
1086      !
1087   END SUBROUTINE rhg_evp_rst
1088
1089   
1090#else
1091   !!----------------------------------------------------------------------
1092   !!   Default option         Empty module           NO SI3 sea-ice model
1093   !!----------------------------------------------------------------------
1094#endif
1095
1096   !!==============================================================================
1097END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.