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

source: NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 @ 15385

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

cleaning some ice routines. No change in sette

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