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/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/ICE/icedyn_rhg_evp.F90 @ 15349

Last change on this file since 15349 was 15349, checked in by cetlod, 2 years ago

dev_PISCO : update with trunk@15348

  • Property svn:keywords set to Id
File size: 65.9 KB
Line 
1MODULE icedyn_rhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_evp  ***
4   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model
8   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
9   !!            3.3  !  2009-05  (G.Garric)    addition of the evp case
10   !!            3.4  !  2011-01  (A. Porter)   dynamical allocation
11   !!            3.5  !  2012-08  (R. Benshila) AGRIF
12   !!            3.6  !  2016-06  (C. Rousset)  Rewriting + landfast ice + mEVP (Bouillon 2013)
13   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
14   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
15   !!----------------------------------------------------------------------
16#if defined key_si3
17   !!----------------------------------------------------------------------
18   !!   'key_si3'                                       SI3 sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
21   !!   rhg_evp_rst     : read/write EVP fields in ice restart
22   !!----------------------------------------------------------------------
23   USE phycst         ! Physical constant
24   USE dom_oce        ! Ocean domain
25   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
26   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
27   USE ice            ! sea-ice: ice variables
28   USE icevar         ! ice_var_sshdyn
29   USE icedyn_rdgrft  ! sea-ice: ice strength
30   USE bdy_oce , ONLY : ln_bdy
31   USE bdyice
32#if defined key_agrif
33   USE agrif_ice_interp
34#endif
35   !
36   USE in_out_manager ! I/O manager
37   USE iom            ! I/O manager library
38   USE lib_mpp        ! MPP library
39   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
40   USE lbclnk         ! lateral boundary conditions (or mpp links)
41   USE prtctl         ! Print control
42
43   USE netcdf         ! NetCDF library for convergence test
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
48   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
49
50   !! for convergence tests
51   INTEGER ::   ncvgid   ! netcdf file id
52   INTEGER ::   nvarid   ! netcdf variable id
53   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   fimask   ! mask at F points for the ice
54   
55   !! * Substitutions
56#  include "do_loop_substitute.h90"
57#  include "domzgr_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
60   !! $Id$
61   !! Software governed by the CeCILL license (see ./LICENSE)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
66      !!-------------------------------------------------------------------
67      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
68      !!                             EVP-C-grid
69      !!
70      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
71      !!  stress and sea-surface slope. Ice-ice interaction is described by
72      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
73      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).
74      !!
75      !!  The points in the C-grid look like this, dear reader
76      !!
77      !!                              (ji,jj)
78      !!                                 |
79      !!                                 |
80      !!                      (ji-1,jj)  |  (ji,jj)
81      !!                             ---------
82      !!                            |         |
83      !!                            | (ji,jj) |------(ji,jj)
84      !!                            |         |
85      !!                             ---------
86      !!                     (ji-1,jj-1)     (ji,jj-1)
87      !!
88      !! ** Inputs  : - wind forcing (stress), oceanic currents
89      !!                ice total volume (vt_i) per unit area
90      !!                snow total volume (vt_s) per unit area
91      !!
92      !! ** Action  : - compute u_ice, v_ice : the components of the
93      !!                sea-ice velocity vector
94      !!              - compute delta_i, shear_i, divu_i, which are inputs
95      !!                of the ice thickness distribution
96      !!
97      !! ** Steps   : 0) compute mask at F point
98      !!              1) Compute ice snow mass, ice strength
99      !!              2) Compute wind, oceanic stresses, mass terms and
100      !!                 coriolis terms of the momentum equation
101      !!              3) Solve the momentum equation (iterative procedure)
102      !!              4) Recompute delta, shear and divergence
103      !!                 (which are inputs of the ITD) & store stress
104      !!                 for the next time step
105      !!              5) Diagnostics including charge ellipse
106      !!
107      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
108      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
109      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
110      !!              (i.e. more stable and better convergence)
111      !!
112      !! References : Hunke and Dukowicz, JPO97
113      !!              Bouillon et al., Ocean Modelling 2009
114      !!              Bouillon et al., Ocean Modelling 2013
115      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
116      !!-------------------------------------------------------------------
117      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
118      INTEGER                 , INTENT(in   ) ::   Kmm                                   ! ocean time level index
119      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
120      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
121      !!
122      INTEGER ::   ji, jj       ! dummy loop indices
123      INTEGER ::   jter         ! local integers
124      !
125      REAL(wp) ::   zrhoco                                              ! rho0 * rn_cio
126      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
127      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
128      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
129      REAl(wp) ::   zbetau, zbetav
130      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
131      REAL(wp) ::   zp_delf, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars
132      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
133      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
134      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
135      !
136      REAL(wp) ::   zintb, zintn                                        ! dummy argument
137      REAL(wp) ::   zfac_x, zfac_y
138      REAL(wp) ::   zshear, zdum1, zdum2
139      !
140      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points
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) ::   zten_i                          ! tension
151      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
152      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
153      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
154      !                                                                 !    ice bottom surface if ice is embedded
155      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
156      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
157      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
158      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
159      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
160      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
161      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
162      !
163      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00, zmsk15
164      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
165      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
166
167      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
168      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
169      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
170      !! --- check convergence
171      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
172      !! --- diags
173      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
175      !! --- SIMIP diags
176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
179      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
180      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
181      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)
182      !! -- advect fields at the rheology time step for the calculation of strength
183      LOGICAL  ::   ll_advups = .FALSE.
184      REAL(wp) ::   zdt_ups
185      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   ztmp
186      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   za_i_ups, zv_i_ups   ! tracers advected upstream
187      !!-------------------------------------------------------------------
188
189      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
190      !
191      ! for diagnostics and convergence tests
192      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
193         zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
194      END_2D
195      IF( nn_rhg_chkcvg > 0 ) THEN
196         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
197            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
198         END_2D
199      ENDIF
200      !
201      !------------------------------------------------------------------------------!
202      ! 0) mask at F points for the ice
203      !------------------------------------------------------------------------------!
204      IF( kt == nit000 ) THEN
205         ! ocean/land mask
206         ALLOCATE( fimask(jpi,jpj) )
207         IF( rn_ishlat == 0._wp ) THEN
208            DO_2D( 0, 0, 0, 0 )
209               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
210            END_2D
211         ELSE
212            DO_2D( 0, 0, 0, 0 )
213               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
214               ! Lateral boundary conditions on velocity (modify fimask)
215               IF( fimask(ji,jj) == 0._wp ) THEN
216                  fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
217                     &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
218               ENDIF
219            END_2D
220         ENDIF
221         CALL lbc_lnk( 'icedyn_rhg_evp', fimask, 'F', 1._wp )
222      ENDIF
223      !------------------------------------------------------------------------------!
224      ! 1) define some variables and initialize arrays
225      !------------------------------------------------------------------------------!
226      zrhoco = rho0 * rn_cio
227
228      ! ecc2: square of yield ellipse eccenticrity
229      ecc2    = rn_ecc * rn_ecc
230      z1_ecc2 = 1._wp / ecc2
231
232      ! alpha parameters (Bouillon 2009)
233      IF( .NOT. ln_aEVP ) THEN
234         zdtevp   = rDt_ice / REAL( nn_nevp )
235         zalph1 =   2._wp * rn_relast * REAL( nn_nevp )
236         zalph2 = zalph1 * z1_ecc2
237
238         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
239         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
240      ELSE
241         zdtevp   = rDt_ice
242         ! zalpha parameters set later on adaptatively
243      ENDIF
244      z1_dtevp = 1._wp / zdtevp
245
246      ! Initialise stress tensor
247      zs1 (:,:) = pstress1_i (:,:)
248      zs2 (:,:) = pstress2_i (:,:)
249      zs12(:,:) = pstress12_i(:,:)
250
251      ! Ice strength
252      CALL ice_strength
253
254      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
255      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
256      ELSE                         ;   zkt = 0._wp
257      ENDIF
258      !
259      !------------------------------------------------------------------------------!
260      ! 2) Wind / ocean stress, mass terms, coriolis terms
261      !------------------------------------------------------------------------------!
262      ! sea surface height
263      !    embedded sea ice: compute representative ice top surface
264      !    non-embedded sea ice: use ocean surface for slope calculation
265      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
266
267      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
268         zm1          = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )  ! Ice/snow mass at U-V points
269         zmf  (ji,jj) = zm1 * ff_t(ji,jj)                            ! Coriolis at T points (m*f)
270         zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin )                   ! dt/m at T points (for alpha and beta coefficients)
271      END_2D
272     
273      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
274
275         ! ice fraction at U-V points
276         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)
277         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)
278
279         ! Ice/snow mass at U-V points
280         zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
281         zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
282         zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
283         zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
284         zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
285
286         ! Ocean currents at U-V points
287         ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
288         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)
289         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)
290
291         ! m/dt
292         zmU_t(ji,jj)    = zmassU * z1_dtevp
293         zmV_t(ji,jj)    = zmassV * z1_dtevp
294
295         ! Drag ice-atm.
296         ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
297         ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
298
299         ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
300         zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
301         zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
302
303         ! masks
304         zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
305         zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
306
307         ! switches
308         IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
309         ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
310         IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
311         ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
312
313      END_2D
314      !
315      !                                  !== Landfast ice parameterization ==!
316      !
317      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
318         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
319            ! ice thickness at U-V points
320            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)
321            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)
322            ! ice-bottom stress at U points
323            zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
324            ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
325            ! ice-bottom stress at V points
326            zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
327            ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
328            ! ice_bottom stress at T points
329            zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) )    ! if grounded icebergs are read: ocean depth = 0
330            tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
331         END_2D
332         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp )
333         !
334      ELSE                               !-- no landfast
335         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
336            ztaux_base(ji,jj) = 0._wp
337            ztauy_base(ji,jj) = 0._wp
338         END_2D
339      ENDIF
340
341      !------------------------------------------------------------------------------!
342      ! 3) Solution of the momentum equation, iterative procedure
343      !------------------------------------------------------------------------------!
344      !
345      !                                               ! ==================== !
346      DO jter = 1 , nn_nevp                           !    loop over jter    !
347         !                                            ! ==================== !
348         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
349         !
350         ! convergence test
351         IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) THEN
352            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
353               zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
354               zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
355            END_2D
356         ENDIF
357
358         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
359         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
360
361            ! shear at F points
362            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)   &
363               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
364               &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
365
366         END_2D
367
368         DO_2D( 0, 0, 0, 0 )
369
370            ! shear**2 at T points (doc eq. A16)
371            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
372               &   + 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)  &
373               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
374
375            ! divergence at T points
376            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
377               &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
378               &    ) * r1_e1e2t(ji,jj)
379            zdiv2 = zdiv * zdiv
380
381            ! tension at T points
382            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)   &
383               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
384               &   ) * r1_e1e2t(ji,jj)
385            zdt2 = zdt * zdt
386
387            ! delta at T points
388            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
389
390            ! P/delta at T points
391            zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
392
393         END_2D
394         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp )
395
396         !
397         DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
398
399            ! divergence at T points (duplication to avoid communications)
400            ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
401            zdiv  = ( (e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj))   &
402               &    + (e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1))   &
403               &    ) * r1_e1e2t(ji,jj)
404
405            ! tension at T points (duplication to avoid communications)
406            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)   &
407               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
408               &   ) * r1_e1e2t(ji,jj)
409
410            ! alpha for aEVP
411            !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
412            !   alpha = beta = sqrt(4*gamma)
413            IF( ln_aEVP ) THEN
414               zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
415               z1_alph1 = 1._wp / ( zalph1 + 1._wp )
416               zalph2   = zalph1
417               z1_alph2 = z1_alph1
418               ! explicit:
419               ! z1_alph1 = 1._wp / zalph1
420               ! z1_alph2 = 1._wp / zalph1
421               ! zalph1 = zalph1 - 1._wp
422               ! zalph2 = zalph1
423            ENDIF
424
425            ! stress at T points (zkt/=0 if landfast)
426            zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1
427            zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
428
429         END_2D
430
431         ! Save beta at T-points for further computations
432         IF( ln_aEVP ) THEN
433            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
434               zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
435            END_2D
436         ENDIF
437
438         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
439
440            ! alpha for aEVP
441            IF( ln_aEVP ) THEN
442               zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) )
443               z1_alph2 = 1._wp / ( zalph2 + 1._wp )
444               ! explicit:
445               ! z1_alph2 = 1._wp / zalph2
446               ! zalph2 = zalph2 - 1._wp
447            ENDIF
448
449            ! P/delta at F points
450            ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
451            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)) )
452
453            ! stress at F points (zkt/=0 if landfast)
454            zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
455
456         END_2D
457
458         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
459         ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
460         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
461            !                   !--- U points
462            zfU(ji,jj) = 0.5_wp * ( (( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
463               &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
464               &                    ) * r1_e2u(ji,jj))                                                                      &
465               &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
466               &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
467               &                  ) * r1_e1e2u(ji,jj)
468            !
469            !                !--- V points
470            zfV(ji,jj) = 0.5_wp * ( (( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
471               &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
472               &                    ) * r1_e1v(ji,jj))                                                                      &
473               &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
474               &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
475               &                  ) * r1_e1e2v(ji,jj)
476            !
477            !                !--- ice currents at U-V point
478            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)
479            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)
480            !
481         END_2D
482         !
483         ! --- Computation of ice velocity --- !
484         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
485         !  Bouillon et al. 2009 (eq 34-35) => stable
486         IF( MOD(jter,2) == 0 ) THEN ! even iterations
487            !
488            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
489               !                 !--- tau_io/(v_oce - v_ice)
490               zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
491                  &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
492               !                 !--- Ocean-to-Ice stress
493               ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
494               !
495               !                 !--- tau_bottom/v_ice
496               zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
497               zTauB = ztauy_base(ji,jj) / zvel
498               !                 !--- OceanBottom-to-Ice stress
499               ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
500               !
501               !                 !--- Coriolis at V-points (energy conserving formulation)
502               zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
503                  &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
504                  &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
505               !
506               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
507               zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
508               !
509               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
510               !                                         1 = sliding friction : TauB < RHS
511               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
512               !
513               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
514                  zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
515                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
516                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
517                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
518                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
519                     &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
520                     &                                    ) / ( zbetav + 1._wp )                                              &
521                     &             ) * 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
522                     &           )   * zmsk00y(ji,jj)
523               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
524                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
525                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
526                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
527                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
528                     &             ) * 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
529                     &            )  * zmsk00y(ji,jj)
530               ENDIF
531            END_2D
532            IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
533            !
534            DO_2D( 0, 0, 0, 0 )
535               !                 !--- tau_io/(u_oce - u_ice)
536               zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
537                  &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
538               !                 !--- Ocean-to-Ice stress
539               ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
540               !
541               !                 !--- tau_bottom/u_ice
542               zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
543               zTauB = ztaux_base(ji,jj) / zvel
544               !                 !--- OceanBottom-to-Ice stress
545               ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
546               !
547               !                 !--- Coriolis at U-points (energy conserving formulation)
548               zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
549                  &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
550                  &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
551               !
552               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
553               zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
554               !
555               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
556               !                                         1 = sliding friction : TauB < RHS
557               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
558               !
559               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
560                  zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
561                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
562                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
563                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
564                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
565                     &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
566                     &                                    ) / ( zbetau + 1._wp )                                              &
567                     &             ) * 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
568                     &           )   * zmsk00x(ji,jj)
569               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
570                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
571                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
572                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
573                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
574                     &             ) * 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
575                     &           )   * zmsk00x(ji,jj)
576               ENDIF
577            END_2D
578            IF( nn_hls == 1 ) THEN   ;   CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
579            ELSE                     ;   CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
580            ENDIF
581            !
582         ELSE ! odd iterations
583            !
584            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
585               !                 !--- tau_io/(u_oce - u_ice)
586               zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
587                  &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
588               !                 !--- Ocean-to-Ice stress
589               ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
590               !
591               !                 !--- tau_bottom/u_ice
592               zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
593               zTauB = ztaux_base(ji,jj) / zvel
594               !                 !--- OceanBottom-to-Ice stress
595               ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
596               !
597               !                 !--- Coriolis at U-points (energy conserving formulation)
598               zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
599                  &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
600                  &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
601               !
602               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
603               zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
604               !
605               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
606               !                                         1 = sliding friction : TauB < RHS
607               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
608               !
609               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
610                  zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
611                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
612                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
613                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
614                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
615                     &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
616                     &                                    ) / ( zbetau + 1._wp )                                              &
617                     &             ) * 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
618                     &           )   * zmsk00x(ji,jj)
619               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
620                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity
621                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
622                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
623                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
624                     &             ) * 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
625                     &           )   * zmsk00x(ji,jj)
626               ENDIF
627            END_2D
628            IF( nn_hls == 1 )   CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp )
629            !
630            DO_2D( 0, 0, 0, 0 )
631               !                 !--- tau_io/(v_oce - v_ice)
632               zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
633                  &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
634               !                 !--- Ocean-to-Ice stress
635               ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
636               !
637               !                 !--- tau_bottom/v_ice
638               zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
639               zTauB = ztauy_base(ji,jj) / zvel
640               !                 !--- OceanBottom-to-Ice stress
641               ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
642               !
643               !                 !--- Coriolis at v-points (energy conserving formulation)
644               zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
645                  &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
646                  &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
647               !
648               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
649               zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
650               !
651               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
652               !                                         1 = sliding friction : TauB < RHS
653               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
654               !
655               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
656                  zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
657                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
658                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
659                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
660                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
661                     &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
662                     &                                    ) / ( zbetav + 1._wp )                                              &
663                     &             ) * 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
664                     &           )   * zmsk00y(ji,jj)
665               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
666                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity
667                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
668                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
669                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
670                     &             ) * 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
671                     &           )   * zmsk00y(ji,jj)
672               ENDIF
673            END_2D
674            IF( nn_hls == 1 ) THEN   ;   CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp )
675            ELSE                     ;   CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp )
676            ENDIF
677            !
678         ENDIF
679         !
680#if defined key_agrif
681!!       CALL agrif_interp_ice( 'U', jter, nn_nevp )
682!!       CALL agrif_interp_ice( 'V', jter, nn_nevp )
683         CALL agrif_interp_ice( 'U' )
684         CALL agrif_interp_ice( 'V' )
685#endif
686         IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
687         IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
688         !
689         ! convergence test
690         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
691         !
692         !
693         ! --- change strength according to advected a_i and v_i (upstream for now) --- !
694         IF( ll_advups .AND. ln_str_H79 ) THEN
695            !
696            IF( jter == 1 ) THEN                               ! init
697               ALLOCATE( za_i_ups(jpi,jpj,jpl), zv_i_ups(jpi,jpj,jpl) )
698               ALLOCATE( ztmp(jpi,jpj) )
699               zdt_ups = rDt_ice / REAL( nn_nevp )
700               za_i_ups(:,:,:) = a_i(:,:,:)
701               zv_i_ups(:,:,:) = v_i(:,:,:)
702            ELSE
703               CALL lbc_lnk( 'icedyn_rhg_evp', za_i_ups, 'T', 1.0_wp, zv_i_ups, 'T', 1.0_wp )               
704            ENDIF
705            !
706            CALL rhg_upstream( zdt_ups, u_ice, v_ice, za_i_ups )   ! upstream advection: a_i
707            CALL rhg_upstream( zdt_ups, u_ice, v_ice, zv_i_ups )   ! upstream advection: v_i
708            !
709            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    ! strength
710               strength(ji,jj) = rn_pstar * SUM( zv_i_ups(ji,jj,:) ) * EXP( -rn_crhg * ( 1._wp - SUM( za_i_ups(ji,jj,:) ) ) )
711            END_2D
712            IF( nn_hls == 1 )  CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp )
713            !
714            DO_2D( 0, 0, 0, 0 )                                ! strength smoothing
715               IF( SUM( za_i_ups(ji,jj,:) ) > 0._wp ) THEN
716                  ztmp(ji,jj) = ( 4._wp * strength(ji,jj) + strength(ji-1,jj  ) + strength(ji+1,jj  ) &
717                     &                                    + strength(ji  ,jj-1) + strength(ji  ,jj+1) &
718                     &          ) / ( 4._wp + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
719               ELSE
720                  ztmp(ji,jj) = 0._wp
721               ENDIF
722            END_2D
723            DO_2D( 0, 0, 0, 0 )
724               strength(ji,jj) = ztmp(ji,jj)
725            END_2D
726            !
727            IF( jter == nn_nevp ) THEN
728               DEALLOCATE( za_i_ups, zv_i_ups )
729               DEALLOCATE( ztmp )
730            ENDIF
731         ENDIF
732         !                                                ! ==================== !
733      END DO                                              !  end loop over jter  !
734      !                                                   ! ==================== !
735      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta )
736      !
737      IF( ll_advups .AND. ln_str_H79 )   CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp )
738      !
739      !------------------------------------------------------------------------------!
740      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
741      !------------------------------------------------------------------------------!
742      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
743
744         ! shear at F points
745         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)   &
746            &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
747            &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
748
749      END_2D
750
751      DO_2D( 0, 0, 0, 0 )   ! no vector loop
752
753         ! tension**2 at T points
754         zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
755            &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
756            &   ) * r1_e1e2t(ji,jj)
757         zdt2 = zdt * zdt
758
759         zten_i(ji,jj) = zdt
760
761         ! shear**2 at T points (doc eq. A16)
762         zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
763            &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
764            &   ) * 0.25_wp * r1_e1e2t(ji,jj)
765
766         ! shear at T points
767         pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
768
769         ! divergence at T points
770         pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
771            &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
772            &             ) * r1_e1e2t(ji,jj)
773
774         ! delta at T points
775         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
776         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
777         pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
778
779      END_2D
780      CALL lbc_lnk( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, &
781         &                            zs1     , 'T', 1._wp, zs2    , 'T', 1._wp, zs12    , 'F', 1._wp )
782
783      ! --- Store the stress tensor for the next time step --- !
784      pstress1_i (:,:) = zs1 (:,:)
785      pstress2_i (:,:) = zs2 (:,:)
786      pstress12_i(:,:) = zs12(:,:)
787      !
788
789      !------------------------------------------------------------------------------!
790      ! 5) diagnostics
791      !------------------------------------------------------------------------------!
792      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
793      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
794         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
795         !
796         CALL lbc_lnk( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, &
797            &                            ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, &
798            &                            ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
799         !
800         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
801         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
802         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
803         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
804         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
805         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
806      ENDIF
807
808      ! --- divergence, shear and strength --- !
809      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
810      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
811      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
812      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
813
814      ! --- Stress tensor invariants (SIMIP diags) --- !
815      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
816         !
817         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
818         !
819         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
820
821            ! Ice stresses
822            ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
823            ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
824            ! I know, this can be confusing...
825            zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )
826            zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
827            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
828            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
829
830            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
831            zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure
832            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
833
834         END_2D
835         !
836         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
837         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
838         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
839
840         DEALLOCATE ( zsig_I, zsig_II )
841
842      ENDIF
843
844      ! --- Normalized stress tensor principal components --- !
845      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
846      ! Recommendation 1 : we use ice strength, not replacement pressure
847      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
848      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
849         !
850         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
851         !
852         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
853
854            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
855            !                        and **deformations** at current iterates
856            !                        following Lemieux & Dupont (2020)
857            zfac             =   zp_delt(ji,jj)
858            zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
859            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
860            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
861
862            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
863            zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                            ! 1st stress invariant, aka average normal stress, aka negative pressure
864            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )   ! 2nd  ''       '', aka maximum shear stress
865
866            ! Normalized  principal stresses (used to display the ellipse)
867            z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
868            zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
869            zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
870         END_2D
871         !
872         CALL iom_put( 'sig1_pnorm' , zsig1_p )
873         CALL iom_put( 'sig2_pnorm' , zsig2_p )
874
875         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
876
877      ENDIF
878
879      ! --- SIMIP --- !
880      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
881         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
882         !
883         CALL lbc_lnk( 'icedyn_rhg_evp', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &
884            &                            zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp )
885
886         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
887         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
888         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
889         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
890         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
891         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
892      ENDIF
893
894      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
895         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
896         !
897         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
898            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
899         !
900         DO_2D( 0, 0, 0, 0 )
901            ! 2D ice mass, snow mass, area transport arrays (X, Y)
902            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
903            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
904
905            zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
906            zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
907
908            zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
909            zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
910
911            zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
912            zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
913
914         END_2D
915
916         CALL lbc_lnk( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &
917            &                            zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, &
918            &                            zdiag_xatrp    , 'U', -1.0_wp, zdiag_yatrp    , 'V', -1.0_wp )
919
920         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
921         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport
922         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
923         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
924         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
925         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
926
927         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
928            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
929
930      ENDIF
931      !
932      ! --- convergence tests --- !
933      IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
934         IF( iom_use('uice_cvg') ) THEN
935            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
936               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
937                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
938            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
939               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
940                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
941            ENDIF
942         ENDIF
943      ENDIF
944      !
945   END SUBROUTINE ice_dyn_rhg_evp
946
947
948   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 )
949      !!----------------------------------------------------------------------
950      !!                    ***  ROUTINE rhg_cvg  ***
951      !!
952      !! ** Purpose :   check convergence of oce rheology
953      !!
954      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity
955      !!                during the sub timestepping of rheology so as:
956      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) )
957      !!                This routine is called every sub-iteration, so it is cpu expensive
958      !!
959      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)
960      !!----------------------------------------------------------------------
961      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
962      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
963      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmsk15
964      !!
965      INTEGER           ::   it, idtime, istatus
966      INTEGER           ::   ji, jj          ! dummy loop indices
967      REAL(wp)          ::   zresm           ! local real
968      CHARACTER(len=20) ::   clname
969      !!----------------------------------------------------------------------
970
971      ! create file
972      IF( kt == nit000 .AND. kiter == 1 ) THEN
973         !
974         IF( lwp ) THEN
975            WRITE(numout,*)
976            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
977            WRITE(numout,*) '~~~~~~~'
978         ENDIF
979         !
980         IF( lwm ) THEN
981            clname = 'ice_cvg.nc'
982            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
983            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
984            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
985            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid )
986            istatus = NF90_ENDDEF(ncvgid)
987         ENDIF
988         !
989      ENDIF
990
991      ! time
992      it = ( kt - 1 ) * kitermax + kiter
993
994      ! convergence
995      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
996         zresm = 0._wp
997      ELSE
998         zresm = 0._wp
999         DO_2D( 0, 0, 0, 0 )
1000            zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
1001               &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
1002         END_2D
1003         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
1004      ENDIF
1005
1006      IF( lwm ) THEN
1007         ! write variables
1008         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
1009         ! close file
1010         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid)
1011      ENDIF
1012
1013   END SUBROUTINE rhg_cvg
1014
1015
1016   SUBROUTINE rhg_evp_rst( cdrw, kt )
1017      !!---------------------------------------------------------------------
1018      !!                   ***  ROUTINE rhg_evp_rst  ***
1019      !!
1020      !! ** Purpose :   Read or write RHG file in restart file
1021      !!
1022      !! ** Method  :   use of IOM library
1023      !!----------------------------------------------------------------------
1024      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1025      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
1026      !
1027      INTEGER  ::   iter            ! local integer
1028      INTEGER  ::   id1, id2, id3   ! local integers
1029      !!----------------------------------------------------------------------
1030      !
1031      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
1032         !                                   ! ---------------
1033         IF( ln_rstart ) THEN                   !* Read the restart file
1034            !
1035            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
1036            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
1037            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
1038            !
1039            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
1040               CALL iom_get( numrir, jpdom_auto, 'stress1_i' , stress1_i , cd_type = 'T' )
1041               CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' )
1042               CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' )
1043            ELSE                                     ! start rheology from rest
1044               IF(lwp) WRITE(numout,*)
1045               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
1046               stress1_i (:,:) = 0._wp
1047               stress2_i (:,:) = 0._wp
1048               stress12_i(:,:) = 0._wp
1049            ENDIF
1050         ELSE                                   !* Start from rest
1051            IF(lwp) WRITE(numout,*)
1052            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
1053            stress1_i (:,:) = 0._wp
1054            stress2_i (:,:) = 0._wp
1055            stress12_i(:,:) = 0._wp
1056         ENDIF
1057         !
1058      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1059         !                                   ! -------------------
1060         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
1061         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
1062         !
1063         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i )
1064         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i )
1065         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
1066         !
1067      ENDIF
1068      !
1069   END SUBROUTINE rhg_evp_rst
1070
1071   SUBROUTINE rhg_upstream( pdt, pu, pv, pt )
1072      !!---------------------------------------------------------------------
1073      !!                    ***  ROUTINE rhg_upstream  ***
1074      !!
1075      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer
1076      !!----------------------------------------------------------------------
1077      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step
1078      REAL(wp), DIMENSION(:,:  ) , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
1079      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   pt               ! tracer fields
1080      !
1081      INTEGER  ::   ji, jj, jl    ! dummy loop indices
1082      REAL(wp) ::   ztra          ! local scalar
1083      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups   ! upstream fluxes
1084      !!----------------------------------------------------------------------
1085      DO jl = 1, jpl
1086         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
1087            zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + &
1088               &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj  ,jl)
1089            zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + &
1090               &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj+1,jl)
1091         END_2D
1092         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
1093            ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )
1094            !
1095            pt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
1096         END_2D
1097      END DO
1098      !
1099   END SUBROUTINE rhg_upstream
1100
1101#else
1102   !!----------------------------------------------------------------------
1103   !!   Default option         Empty module           NO SI3 sea-ice model
1104   !!----------------------------------------------------------------------
1105#endif
1106
1107   !!==============================================================================
1108END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.