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

Last change on this file since 13549 was 13549, checked in by clem, 5 months ago

4.0-HEAD: remove one communication in rheology, so about 100 overall.

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