New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rhg_evp.F90 in NEMO/branches/2020/SI3-03_VP_rheology/src/ICE – NEMO

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

Last change on this file since 15516 was 15516, checked in by vancop, 8 months ago

fixes in diagnostics

  • Property svn:keywords set to Id
File size: 64.4 KB
Line 
1MODULE icedyn_rhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_evp  ***
4   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model
8   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
9   !!            3.3  !  2009-05  (G.Garric)    addition of the evp case
10   !!            3.4  !  2011-01  (A. Porter)   dynamical allocation
11   !!            3.5  !  2012-08  (R. Benshila) AGRIF
12   !!            3.6  !  2016-06  (C. Rousset)  Rewriting + landfast ice + mEVP (Bouillon 2013)
13   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
14   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
15   !!----------------------------------------------------------------------
16#if defined key_si3
17   !!----------------------------------------------------------------------
18   !!   'key_si3'                                       SI3 sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
21   !!   rhg_evp_rst     : read/write EVP fields in ice restart
22   !!----------------------------------------------------------------------
23   USE phycst         ! Physical constant
24   USE dom_oce        ! Ocean domain
25   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
26   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
27   USE ice            ! sea-ice: ice variables
28   USE icevar         ! ice_var_sshdyn
29   USE icedyn_rdgrft  ! sea-ice: ice strength
30   USE bdy_oce , ONLY : ln_bdy 
31   USE bdyice 
32#if defined key_agrif
33   USE agrif_ice_interp
34#endif
35   !
36   USE in_out_manager ! I/O manager
37   USE iom            ! I/O manager library
38   USE lib_mpp        ! MPP library
39   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
40   USE lbclnk         ! lateral boundary conditions (or mpp links)
41   USE prtctl         ! Print control
42
43   USE netcdf         ! NetCDF library for convergence test
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
48   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
49
50   !! * Substitutions
51#  include "vectopt_loop_substitute.h90"
52
53   !! for convergence tests
54   INTEGER ::   ncvgid   ! netcdf file id
55   INTEGER ::   nvarid   ! netcdf variable id
56   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
57   !!----------------------------------------------------------------------
58   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
59   !! $Id$
60   !! Software governed by the CeCILL license (see ./LICENSE)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
65      !!-------------------------------------------------------------------
66      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
67      !!                             EVP-C-grid
68      !!
69      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
70      !!  stress and sea-surface slope. Ice-ice interaction is described by
71      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
72      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
73      !!
74      !!  The points in the C-grid look like this, dear reader
75      !!
76      !!                              (ji,jj)
77      !!                                 |
78      !!                                 |
79      !!                      (ji-1,jj)  |  (ji,jj)
80      !!                             ---------   
81      !!                            |         |
82      !!                            | (ji,jj) |------(ji,jj)
83      !!                            |         |
84      !!                             ---------   
85      !!                     (ji-1,jj-1)     (ji,jj-1)
86      !!
87      !! ** Inputs  : - wind forcing (stress), oceanic currents
88      !!                ice total volume (vt_i) per unit area
89      !!                snow total volume (vt_s) per unit area
90      !!
91      !! ** Action  : - compute u_ice, v_ice : the components of the
92      !!                sea-ice velocity vector
93      !!              - compute delta_i, shear_i, divu_i, which are inputs
94      !!                of the ice thickness distribution
95      !!
96      !! ** Steps   : 0) compute mask at F point
97      !!              1) Compute ice snow mass, ice strength
98      !!              2) Compute wind, oceanic stresses, mass terms and
99      !!                 coriolis terms of the momentum equation
100      !!              3) Solve the momentum equation (iterative procedure)
101      !!              4) Recompute delta, shear and divergence
102      !!                 (which are inputs of the ITD) & store stress
103      !!                 for the next time step
104      !!              5) Diagnostics including charge ellipse
105      !!
106      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
107      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
108      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
109      !!              (i.e. more stable and better convergence)
110      !!
111      !! References : Hunke and Dukowicz, JPO97
112      !!              Bouillon et al., Ocean Modelling 2009
113      !!              Bouillon et al., Ocean Modelling 2013
114      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
115      !!-------------------------------------------------------------------
116      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
117      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
118      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
119      !!
120      INTEGER ::   ji, jj       ! dummy loop indices
121      INTEGER ::   jter         ! local integers
122      !
123      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
124      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
125      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
126      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
127      REAl(wp) ::   zbetau, zbetav
128      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
129      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
130      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
131      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
132      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
133      !
134      REAL(wp) ::   zintb, zintn                                        ! dummy argument
135      REAL(wp) ::   zfac_x, zfac_y
136      REAL(wp) ::   zdum1, zdum2
137      !
138      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
139      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltat                         ! Delta at T-points
140      REAL(wp), DIMENSION(jpi,jpj) ::   ztens, zshear                   ! Tension and shear
141      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
142      !
143      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
144      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
145      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
146      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
147      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
148      !
149      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
150      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
151      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
152      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
153      !                                                                 !    ice bottom surface if ice is embedded   
154      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
155      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
156      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
157      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
158      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
159      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
160      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
161      !
162      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
163      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
164      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
165
166      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
167      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
168      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
169      !! --- check convergence
170      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
171      !! --- diags
172      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
174      !! --- SIMIP diags
175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
179      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
180      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
181      !!-------------------------------------------------------------------
182
183      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
184      !
185      ! for diagnostics and convergence tests
186      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
187      DO jj = 1, jpj
188         DO ji = 1, jpi
189            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
190            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
191         END DO
192      END DO
193      !
194      !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
195      !------------------------------------------------------------------------------!
196      ! 0) mask at F points for the ice
197      !------------------------------------------------------------------------------!
198      ! ocean/land mask
199      DO jj = 1, jpjm1
200         DO ji = 1, jpim1      ! NO vector opt.
201            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
202         END DO
203      END DO
204      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
205
206      ! Lateral boundary conditions on velocity (modify zfmask)
207      DO jj = 2, jpjm1
208         DO ji = fs_2, fs_jpim1   ! vector opt.
209            IF( zfmask(ji,jj) == 0._wp ) THEN
210               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
211                  &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
212            ENDIF
213         END DO
214      END DO
215      DO jj = 2, jpjm1
216         IF( zfmask(1,jj) == 0._wp ) THEN
217            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) )
218         ENDIF
219         IF( zfmask(jpi,jj) == 0._wp ) THEN
220            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) )
221        ENDIF
222      END DO
223      DO ji = 2, jpim1
224         IF( zfmask(ji,1) == 0._wp ) THEN
225            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) )
226         ENDIF
227         IF( zfmask(ji,jpj) == 0._wp ) THEN
228            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) )
229         ENDIF
230      END DO
231      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
232
233      !------------------------------------------------------------------------------!
234      ! 1) define some variables and initialize arrays
235      !------------------------------------------------------------------------------!
236      zrhoco = rau0 * rn_cio 
237
238      ! ecc2: square of yield ellipse eccenticrity
239      ecc2    = rn_ecc * rn_ecc
240      z1_ecc2 = 1._wp / ecc2
241
242      ! alpha parameters (Bouillon 2009)
243      IF( .NOT. ln_aEVP ) THEN
244         zdtevp   = rdt_ice / REAL( nn_nevp )
245         zalph1 =   2._wp * rn_relast * REAL( nn_nevp )
246         zalph2 = zalph1 * z1_ecc2
247
248         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
249         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
250      ELSE
251         zdtevp   = rdt_ice
252         ! zalpha parameters set later on adaptatively
253      ENDIF
254      z1_dtevp = 1._wp / zdtevp
255         
256      ! Initialise stress tensor
257      zs1 (:,:) = pstress1_i (:,:) 
258      zs2 (:,:) = pstress2_i (:,:)
259      zs12(:,:) = pstress12_i(:,:)
260
261      ! Ice strength
262      CALL ice_strength
263
264      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
265      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
266      ELSE                         ;   zkt = 0._wp
267      ENDIF
268      !
269      !------------------------------------------------------------------------------!
270      ! 2) Wind / ocean stress, mass terms, coriolis terms
271      !------------------------------------------------------------------------------!
272      ! sea surface height
273      !    embedded sea ice: compute representative ice top surface
274      !    non-embedded sea ice: use ocean surface for slope calculation
275      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
276
277      DO jj = 2, jpjm1
278         DO ji = fs_2, fs_jpim1
279
280            ! ice fraction at U-V points
281            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)
282            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)
283
284            ! Ice/snow mass at U-V points
285            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
286            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
287            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
288            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
289            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
290
291            ! Ocean currents at U-V points
292            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)
293            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)
294
295            ! Coriolis at T points (m*f)
296            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
297
298            ! dt/m at T points (for alpha and beta coefficients)
299            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
300           
301            ! m/dt
302            zmU_t(ji,jj)    = zmassU * z1_dtevp
303            zmV_t(ji,jj)    = zmassV * z1_dtevp
304           
305            ! Drag ice-atm.
306            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
307            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
308
309            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
310            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
311            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
312
313            ! masks
314            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
315            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
316
317            ! switches
318            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
319            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
320            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
321            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
322
323         END DO
324      END DO
325      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. )
326      !
327      !                                  !== Landfast ice parameterization ==!
328      !
329      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
330         DO jj = 2, jpjm1
331            DO ji = fs_2, fs_jpim1
332               ! ice thickness at U-V points
333               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)
334               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)
335               ! ice-bottom stress at U points
336               zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj)
337               ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
338               ! ice-bottom stress at V points
339               zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj)
340               ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
341               ! ice_bottom stress at T points
342               zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj)
343               tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
344            END DO
345         END DO
346         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. )
347         !
348      ELSE                               !-- no landfast
349         DO jj = 2, jpjm1
350            DO ji = fs_2, fs_jpim1
351               ztaux_base(ji,jj) = 0._wp
352               ztauy_base(ji,jj) = 0._wp
353            END DO
354         END DO
355      ENDIF
356
357      !------------------------------------------------------------------------------!
358      ! 3) Solution of the momentum equation, iterative procedure
359      !------------------------------------------------------------------------------!
360      !
361      !                                               ! ==================== !
362      DO jter = 1 , nn_nevp                           !    loop over jter    !
363         !                                            ! ==================== !       
364         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
365         !
366         ! convergence test
367         IF( ln_rhg_chkcvg ) THEN
368            DO jj = 1, jpj
369               DO ji = 1, jpi
370                  zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
371                  zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
372               END DO
373            END DO
374         ENDIF
375
376         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
377         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
378            DO ji = 1, jpim1
379
380               ! shear at F points
381               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)   &
382                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
383                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
384
385            END DO
386         END DO
387         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
388
389         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
390            DO ji = 2, jpi ! no vector loop
391
392               ! shear**2 at T points (doc eq. A16)
393               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
394                  &   + 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)  &
395                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
396             
397               ! divergence at T points
398               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
399                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
400                  &    ) * r1_e1e2t(ji,jj)
401               zdiv2 = zdiv * zdiv
402               
403               ! tension at T points
404               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)   &
405                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
406                  &   ) * r1_e1e2t(ji,jj)
407               zdt2 = zdt * zdt
408               
409               ! delta at T points
410               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
411               
412               zdeltat(ji,jj) = zdelta  ! store delta at previous iterate (for ellipse computation)
413
414               ! P/delta at T points
415               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
416
417               ! alpha for aEVP
418               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
419               !   alpha = beta = sqrt(4*gamma)
420               IF( ln_aEVP ) THEN
421                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
422                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
423                  zalph2   = zalph1
424                  z1_alph2 = z1_alph1
425                  ! explicit:
426                  ! z1_alph1 = 1._wp / zalph1
427                  ! z1_alph2 = 1._wp / zalph1
428                  ! zalph1 = zalph1 - 1._wp
429                  ! zalph2 = zalph1
430               ENDIF
431               
432               ! stress at T points (zkt/=0 if landfast)
433               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
434               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
435             
436            END DO
437         END DO
438         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
439
440         ! Save beta at T-points for further computations
441         IF( ln_aEVP ) THEN
442            DO jj = 1, jpj
443               DO ji = 1, jpi
444                  zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
445               END DO
446            END DO
447         ENDIF
448         
449         DO jj = 1, jpjm1
450            DO ji = 1, jpim1
451
452               ! alpha for aEVP
453               IF( ln_aEVP ) THEN
454                  zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) )
455                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
456                  ! explicit:
457                  ! z1_alph2 = 1._wp / zalph2
458                  ! zalph2 = zalph2 - 1._wp
459               ENDIF
460               
461               ! P/delta at F points
462               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) )
463               
464               ! stress at F points (zkt/=0 if landfast)
465               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
466
467            END DO
468         END DO
469
470         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
471         DO jj = 2, jpjm1
472            DO ji = fs_2, fs_jpim1               
473               !                   !--- U points
474               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
475                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
476                  &                    ) * r1_e2u(ji,jj)                                                                      &
477                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
478                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
479                  &                  ) * r1_e1e2u(ji,jj)
480               !
481               !                !--- V points
482               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
483                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
484                  &                    ) * r1_e1v(ji,jj)                                                                      &
485                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
486                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
487                  &                  ) * r1_e1e2v(ji,jj)
488               !
489               !                !--- ice currents at U-V point
490               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)
491               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)
492               !
493            END DO
494         END DO
495         !
496         ! --- Computation of ice velocity --- !
497         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
498         !  Bouillon et al. 2009 (eq 34-35) => stable
499         IF( MOD(jter,2) == 0 ) THEN ! even iterations
500            !
501            DO jj = 2, jpjm1
502               DO ji = fs_2, fs_jpim1
503                  !                 !--- tau_io/(v_oce - v_ice)
504                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
505                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
506                  !                 !--- Ocean-to-Ice stress
507                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
508                  !
509                  !                 !--- tau_bottom/v_ice
510                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
511                  zTauB = ztauy_base(ji,jj) / zvel
512                  !                 !--- OceanBottom-to-Ice stress
513                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
514                  !
515                  !                 !--- Coriolis at V-points (energy conserving formulation)
516                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
517                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
518                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
519                  !
520                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
521                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
522                  !
523                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
524                  !                                         1 = sliding friction : TauB < RHS
525                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
526                  !
527                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
528                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
529                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
530                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
531                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
532                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
533                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
534                        &                                    ) / ( zbetav + 1._wp )                                              &
535                        &             ) * 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
536                        &           )   * zmsk00y(ji,jj)
537                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
538                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
539                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
540                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
541                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
542                        &             ) * 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
543                        &            )  * zmsk00y(ji,jj)
544                  ENDIF
545               END DO
546            END DO
547            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
548            !
549#if defined key_agrif
550!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
551            CALL agrif_interp_ice( 'V' )
552#endif
553            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
554            !
555            DO jj = 2, jpjm1
556               DO ji = fs_2, fs_jpim1         
557                  !                 !--- tau_io/(u_oce - u_ice)
558                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
559                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
560                  !                 !--- Ocean-to-Ice stress
561                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
562                  !
563                  !                 !--- tau_bottom/u_ice
564                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
565                  zTauB = ztaux_base(ji,jj) / zvel
566                  !                 !--- OceanBottom-to-Ice stress
567                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
568                  !
569                  !                 !--- Coriolis at U-points (energy conserving formulation)
570                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
571                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
572                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
573                  !
574                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
575                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
576                  !
577                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
578                  !                                         1 = sliding friction : TauB < RHS
579                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
580                  !
581                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
582                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
583                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
584                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
585                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
586                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
587                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
588                        &                                    ) / ( zbetau + 1._wp )                                              &
589                        &             ) * 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
590                        &           )   * zmsk00x(ji,jj)
591                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
592                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
593                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
594                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
595                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
596                        &             ) * 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
597                        &           )   * zmsk00x(ji,jj)
598                  ENDIF
599               END DO
600            END DO
601            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
602            !
603#if defined key_agrif
604!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
605            CALL agrif_interp_ice( 'U' )
606#endif
607            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
608            !
609         ELSE ! odd iterations
610            !
611            DO jj = 2, jpjm1
612               DO ji = fs_2, fs_jpim1
613                  !                 !--- tau_io/(u_oce - u_ice)
614                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
615                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
616                  !                 !--- Ocean-to-Ice stress
617                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
618                  !
619                  !                 !--- tau_bottom/u_ice
620                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
621                  zTauB = ztaux_base(ji,jj) / zvel
622                  !                 !--- OceanBottom-to-Ice stress
623                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
624                  !
625                  !                 !--- Coriolis at U-points (energy conserving formulation)
626                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
627                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
628                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
629                  !
630                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
631                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
632                  !
633                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
634                  !                                         1 = sliding friction : TauB < RHS
635                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
636                  !
637                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
638                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
639                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
640                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
641                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
642                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
643                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
644                        &                                    ) / ( zbetau + 1._wp )                                              &
645                        &             ) * 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
646                        &           )   * zmsk00x(ji,jj)
647                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
648                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
649                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
650                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
651                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
652                        &             ) * 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
653                        &           )   * zmsk00x(ji,jj)
654                  ENDIF
655               END DO
656            END DO
657            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
658            !
659#if defined key_agrif
660!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
661            CALL agrif_interp_ice( 'U' )
662#endif
663            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
664            !
665            DO jj = 2, jpjm1
666               DO ji = fs_2, fs_jpim1
667                  !                 !--- tau_io/(v_oce - v_ice)
668                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
669                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
670                  !                 !--- Ocean-to-Ice stress
671                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
672                  !
673                  !                 !--- tau_bottom/v_ice
674                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
675                  zTauB = ztauy_base(ji,jj) / zvel
676                  !                 !--- OceanBottom-to-Ice stress
677                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
678                  !
679                  !                 !--- Coriolis at v-points (energy conserving formulation)
680                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
681                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
682                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
683                  !
684                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
685                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
686                  !
687                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
688                  !                                         1 = sliding friction : TauB < RHS
689                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
690                  !
691                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
692                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
693                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
694                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
695                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
696                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
697                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
698                        &                                    ) / ( zbetav + 1._wp )                                              & 
699                        &             ) * 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
700                        &           )   * zmsk00y(ji,jj)
701                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
702                                     !---> MV : comment above is misleading as EVP is a purely explicit method
703                                     !---> MV : I found code below quite unreadable :-)
704                                     !---> MV : two maskings + landfast / no landfast options, I think this is a bit too much
705                                     !---> MV : I would suggest to split at least the landfast / masking steps if possible
706                                     !---> MV : Should also comment why 1% of ocean velocity is assumed - not sure it makes sense btw
707                                     !---> MV : I would say 100% of ocean current makes much more sense for free-drift
708                                     !---> MV : of very low concentration sea ice..
709                                     !---> MV : but maybe it is a numerical trick... in brief it's n
710                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
711                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
712                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
713                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
714                        &             ) * 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
715                        &           )   * zmsk00y(ji,jj)
716                  ENDIF
717               END DO
718            END DO
719            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
720            !
721#if defined key_agrif
722!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
723            CALL agrif_interp_ice( 'V' )
724#endif
725            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
726            !
727         ENDIF
728
729         ! convergence test
730         IF( ln_rhg_chkcvg )   CALL rhg_cvg_evp( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
731         !
732         !                                                ! ==================== !
733      END DO                                              !  end loop over jter  !
734      !                                                   ! ==================== !
735      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta )
736      !
737      !------------------------------------------------------------------------------!
738      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
739      !------------------------------------------------------------------------------!
740      DO jj = 1, jpjm1
741         DO ji = 1, jpim1
742
743            ! shear at F points
744            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)   &
745               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
746               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
747
748         END DO
749      END DO           
750     
751      DO jj = 2, jpjm1
752         DO ji = 2, jpim1 ! no vector loop
753           
754            ! tension**2 at T points
755            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)   &
756               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
757               &   ) * r1_e1e2t(ji,jj)
758            zdt2 = zdt * zdt
759           
760            ztens(ji,jj) = zdt
761           
762            ! shear**2 at T points (doc eq. A16)
763            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
764               &   + 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)  &
765               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
766
767            zshear(ji,jj) = SQRT( zds2 )
768           
769            ! maximum shear rate at T points (includes tension, output only)
770            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
771
772            ! divergence at T points
773            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
774               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
775               &             ) * r1_e1e2t(ji,jj)
776           
777            ! delta at T points
778            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
779            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
780            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
781
782         END DO
783      END DO
784      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
785     
786      ! --- Store the stress tensor for the next time step --- !
787      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
788      pstress1_i (:,:) = zs1 (:,:)
789      pstress2_i (:,:) = zs2 (:,:)
790      pstress12_i(:,:) = zs12(:,:)
791      !
792
793      !------------------------------------------------------------------------------!
794      ! 5) diagnostics
795      !------------------------------------------------------------------------------!
796      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
797      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
798         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
799         !
800         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
801            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
802         !
803         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
804         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
805         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
806         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
807         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
808         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
809      ENDIF
810       
811      ! --- divergence, shear and strength --- !
812      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
813      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
814      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
815      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
816
817      ! --- Stress tensor invariants (SIMIP diags) --- !
818      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
819         !
820         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
821         !         
822         DO jj = 2, jpj - 1
823            DO ji = 2, jpi - 1
824           
825               ! Ice stresses
826               ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
827               ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
828               ! I know, this can be confusing...
829               zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 
830               zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
831               zsig2            =   zfac * z1_ecc2 * ztens(ji,jj)
832               zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
833               
834               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
835               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                       ! 1st stress invariant, aka average normal stress, aka negative pressure
836               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )   ! 2nd  ''       '', aka maximum shear stress
837               
838            END DO
839         END DO
840
841         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
842         
843         !
844         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
845         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,   zsig_I(:,:)  * zmsk00(:,:) ) ! Normal stress
846         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' ,   zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
847         
848         DEALLOCATE ( zsig_I, zsig_II )
849         
850      ENDIF
851
852      ! --- Normalized stress tensor principal components --- !
853      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
854      ! Recommendation 1 : we use ice strength, not replacement pressure
855      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
856      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
857         !
858         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
859         !         
860         DO jj = 2, jpj - 1
861            DO ji = 2, jpi - 1
862           
863               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
864               !                        and **deformations** at current iterates
865               !                        following Lemieux & Dupont (2020)
866               zfac             =   zp_delt(ji,jj)
867               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltat(ji,jj) )
868               zsig2            =   zfac * z1_ecc2 * ztens(ji,jj)
869               zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
870               
871               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
872               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                ! 1st stress invariant, aka average normal stress, aka negative pressure
873               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )     ! 2nd  ''       '', aka maximum shear stress
874
875               ! Normalized  principal stresses (used to display the ellipse)
876               z1_strength      =   1._wp / MAX ( 1.0, strength(ji,jj) )
877               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
878               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
879            END DO
880         END DO
881         
882         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
883               
884         !
885         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
886         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
887
888         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
889         
890      ENDIF
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( ln_rhg_chkcvg ) 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_evp( kt, kiter, kitermax, pu, pv, pub, pvb )
965      !!----------------------------------------------------------------------
966      !!                    ***  ROUTINE rhg_cvg_evp  ***
967      !!                     
968      !! ** Purpose :   check convergence of evp ice 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_evp : 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 - nit000 ) * 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 .AND. kiter == kitermax )    istatus = NF90_CLOSE( ncvgid )
1029      ENDIF
1030     
1031   END SUBROUTINE rhg_cvg_evp
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.