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

Last change on this file since 10425 was 10425, checked in by smasson, 20 months ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 57.7 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   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
47   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
53   !! $Id$
54   !! Software governed by the CeCILL license (see ./LICENSE)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
59      !!-------------------------------------------------------------------
60      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
61      !!                             EVP-C-grid
62      !!
63      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
64      !!  stress and sea-surface slope. Ice-ice interaction is described by
65      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
66      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
67      !!
68      !!  The points in the C-grid look like this, dear reader
69      !!
70      !!                              (ji,jj)
71      !!                                 |
72      !!                                 |
73      !!                      (ji-1,jj)  |  (ji,jj)
74      !!                             ---------   
75      !!                            |         |
76      !!                            | (ji,jj) |------(ji,jj)
77      !!                            |         |
78      !!                             ---------   
79      !!                     (ji-1,jj-1)     (ji,jj-1)
80      !!
81      !! ** Inputs  : - wind forcing (stress), oceanic currents
82      !!                ice total volume (vt_i) per unit area
83      !!                snow total volume (vt_s) per unit area
84      !!
85      !! ** Action  : - compute u_ice, v_ice : the components of the
86      !!                sea-ice velocity vector
87      !!              - compute delta_i, shear_i, divu_i, which are inputs
88      !!                of the ice thickness distribution
89      !!
90      !! ** Steps   : 0) compute mask at F point
91      !!              1) Compute ice snow mass, ice strength
92      !!              2) Compute wind, oceanic stresses, mass terms and
93      !!                 coriolis terms of the momentum equation
94      !!              3) Solve the momentum equation (iterative procedure)
95      !!              4) Recompute delta, shear and divergence
96      !!                 (which are inputs of the ITD) & store stress
97      !!                 for the next time step
98      !!              5) Diagnostics including charge ellipse
99      !!
100      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
101      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
102      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
103      !!              (i.e. more stable and better convergence)
104      !!
105      !! References : Hunke and Dukowicz, JPO97
106      !!              Bouillon et al., Ocean Modelling 2009
107      !!              Bouillon et al., Ocean Modelling 2013
108      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
109      !!-------------------------------------------------------------------
110      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
111      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
113      !!
114      INTEGER ::   ji, jj       ! dummy loop indices
115      INTEGER ::   jter         ! local integers
116      !
117      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
118      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
119      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
120      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
121      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
122      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
123      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars
124      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
125      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
126      !
127      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
128      REAL(wp) ::   zintb, zintn                                        ! dummy argument
129      REAL(wp) ::   zfac_x, zfac_y
130      REAL(wp) ::   zshear, zdum1, zdum2
131      !
132      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
133      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
134      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
135      !
136      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
137      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
138      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
139      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
140      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
141      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param)
142      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
143      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
144      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
145      !
146      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
147      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
148!!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
149      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
150      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
151      !                                                                 !    ice bottom surface if ice is embedded   
152      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
153      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
154      !
155      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
156      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
157      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
158
159      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
160      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
161      !! --- diags
162      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
163      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
164      !! --- SIMIP diags
165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress 
175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
179      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
180      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
181      !!-------------------------------------------------------------------
182
183      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
184      !
185!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
186      !------------------------------------------------------------------------------!
187      ! 0) mask at F points for the ice
188      !------------------------------------------------------------------------------!
189      ! ocean/land mask
190      DO jj = 1, jpjm1
191         DO ji = 1, jpim1      ! NO vector opt.
192            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
193         END DO
194      END DO
195      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
196
197      ! Lateral boundary conditions on velocity (modify zfmask)
198      zwf(:,:) = zfmask(:,:)
199      DO jj = 2, jpjm1
200         DO ji = fs_2, fs_jpim1   ! vector opt.
201            IF( zfmask(ji,jj) == 0._wp ) THEN
202               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
203            ENDIF
204         END DO
205      END DO
206      DO jj = 2, jpjm1
207         IF( zfmask(1,jj) == 0._wp ) THEN
208            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
209         ENDIF
210         IF( zfmask(jpi,jj) == 0._wp ) THEN
211            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
212         ENDIF
213      END DO
214      DO ji = 2, jpim1
215         IF( zfmask(ji,1) == 0._wp ) THEN
216            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
217         ENDIF
218         IF( zfmask(ji,jpj) == 0._wp ) THEN
219            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
220         ENDIF
221      END DO
222      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
223
224      !------------------------------------------------------------------------------!
225      ! 1) define some variables and initialize arrays
226      !------------------------------------------------------------------------------!
227      zrhoco = rau0 * rn_cio 
228
229      ! ecc2: square of yield ellipse eccenticrity
230      ecc2    = rn_ecc * rn_ecc
231      z1_ecc2 = 1._wp / ecc2
232
233      ! Time step for subcycling
234      zdtevp   = rdt_ice / REAL( nn_nevp )
235      z1_dtevp = 1._wp / zdtevp
236
237      ! alpha parameters (Bouillon 2009)
238      IF( .NOT. ln_aEVP ) THEN
239         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
240         zalph2 = zalph1 * z1_ecc2
241
242         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
243         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
244      ENDIF
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      ! scale factors
255      DO jj = 2, jpjm1
256         DO ji = fs_2, fs_jpim1
257            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
258            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
259         END DO
260      END DO
261
262      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
263      IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile
264      ELSE                                               ;   zkt = 0._wp
265      ENDIF
266      !
267      !------------------------------------------------------------------------------!
268      ! 2) Wind / ocean stress, mass terms, coriolis terms
269      !------------------------------------------------------------------------------!
270      ! sea surface height
271      !    embedded sea ice: compute representative ice top surface
272      !    non-embedded sea ice: use ocean surface for slope calculation
273      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
274
275      DO jj = 2, jpjm1
276         DO ji = fs_2, fs_jpim1
277
278            ! ice fraction at U-V points
279            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)
280            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)
281
282            ! Ice/snow mass at U-V points
283            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
284            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
285            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
286            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
287            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
288
289            ! Ocean currents at U-V points
290            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
291               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
292           
293            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
294               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
295
296            ! Coriolis at T points (m*f)
297            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
298
299            ! dt/m at T points (for alpha and beta coefficients)
300            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
301           
302            ! m/dt
303            zmU_t(ji,jj)    = zmassU * z1_dtevp
304            zmV_t(ji,jj)    = zmassV * z1_dtevp
305           
306            ! Drag ice-atm.
307            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
308            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
309
310            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
311            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
312            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
313
314            ! masks
315            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
316            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
317
318            ! switches
319            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
320            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
321
322         END DO
323      END DO
324      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. )
325      !
326      !                                  !== Landfast ice parameterization ==!
327      !
328      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
329         DO jj = 2, jpjm1
330            DO ji = fs_2, fs_jpim1
331               ! ice thickness at U-V points
332               zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
333               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
334               ! ice-bottom stress at U points
335               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
336               zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
337               ! ice-bottom stress at V points
338               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
339               zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
340               ! ice_bottom stress at T points
341               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
342               tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
343            END DO
344         END DO
345         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. )
346         !
347      ELSEIF( ln_landfast_home ) THEN          !-- Home made
348         DO jj = 2, jpjm1
349            DO ji = fs_2, fs_jpim1
350               zTauU_ib(ji,jj) = tau_icebfr(ji,jj)
351               zTauV_ib(ji,jj) = tau_icebfr(ji,jj)
352            END DO
353         END DO
354         !
355      ELSE                                     !-- no landfast
356         DO jj = 2, jpjm1
357            DO ji = fs_2, fs_jpim1
358               zTauU_ib(ji,jj) = 0._wp
359               zTauV_ib(ji,jj) = 0._wp
360            END DO
361         END DO
362      ENDIF
363      IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )
364
365      !------------------------------------------------------------------------------!
366      ! 3) Solution of the momentum equation, iterative procedure
367      !------------------------------------------------------------------------------!
368      !
369      !                                               !----------------------!
370      DO jter = 1 , nn_nevp                           !    loop over jter    !
371         !                                            !----------------------!       
372         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
373         !
374!!$         IF(ln_ctl) THEN   ! Convergence test
375!!$            DO jj = 1, jpjm1
376!!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
377!!$               zv_ice(:,jj) = v_ice(:,jj)
378!!$            END DO
379!!$         ENDIF
380
381         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
382         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
383            DO ji = 1, jpim1
384
385               ! shear at F points
386               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)   &
387                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
388                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
389
390            END DO
391         END DO
392         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
393
394         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
395            DO ji = 2, jpi ! no vector loop
396
397               ! shear**2 at T points (doc eq. A16)
398               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
399                  &   + 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)  &
400                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
401             
402               ! divergence at T points
403               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
404                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
405                  &    ) * r1_e1e2t(ji,jj)
406               zdiv2 = zdiv * zdiv
407               
408               ! tension at T points
409               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)   &
410                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
411                  &   ) * r1_e1e2t(ji,jj)
412               zdt2 = zdt * zdt
413               
414               ! delta at T points
415               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
416
417               ! P/delta at T points
418               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
419
420               ! alpha & beta for aEVP
421               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
422               !   alpha = beta = sqrt(4*gamma)
423               IF( ln_aEVP ) THEN
424                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
425                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
426                  zalph2   = zalph1
427                  z1_alph2 = z1_alph1
428               ENDIF
429               
430               ! stress at T points (zkt/=0 if landfast)
431               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
432               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
433             
434            END DO
435         END DO
436         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
437
438         DO jj = 1, jpjm1
439            DO ji = 1, jpim1
440
441               ! alpha & beta for aEVP
442               IF( ln_aEVP ) THEN
443                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
444                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
445                  zbeta(ji,jj) = zalph2
446               ENDIF
447               
448               ! P/delta at F points
449               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) )
450               
451               ! stress at F points (zkt/=0 if landfast)
452               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
453
454            END DO
455         END DO
456
457         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
458         DO jj = 2, jpjm1
459            DO ji = fs_2, fs_jpim1               
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               !                !--- u_ice at V point
477               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
478                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
479               !
480               !                !--- v_ice at U point
481               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
482                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
483            END DO
484         END DO
485         !
486         ! --- Computation of ice velocity --- !
487         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
488         !  Bouillon et al. 2009 (eq 34-35) => stable
489         IF( MOD(jter,2) == 0 ) THEN ! even iterations
490            !
491            DO jj = 2, jpjm1
492               DO ji = fs_2, fs_jpim1
493                  !                 !--- tau_io/(v_oce - v_ice)
494                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
495                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
496                  !                 !--- Ocean-to-Ice stress
497                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
498                  !
499                  !                 !--- tau_bottom/v_ice
500                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
501                  zTauB = - zTauV_ib(ji,jj) / zvel
502                  !
503                  !                 !--- Coriolis at V-points (energy conserving formulation)
504                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
505                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
506                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
507                  !
508                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
509                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
510                  !
511                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
512                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
513                  !
514                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
515                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
516                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
517                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
518                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
519                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
520                     &           ) * zmaskV(ji,jj)
521                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
522                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
523                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
524                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
525                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
526                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
527                     &            ) * zmaskV(ji,jj)
528                  ENDIF
529               END DO
530            END DO
531            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
532            !
533#if defined key_agrif
534!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
535            CALL agrif_interp_ice( 'V' )
536#endif
537            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
538            !
539            DO jj = 2, jpjm1
540               DO ji = fs_2, fs_jpim1         
541                  !                 !--- tau_io/(u_oce - u_ice)
542                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
543                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
544                  !                 !--- Ocean-to-Ice stress
545                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
546                  !
547                  !                 !--- tau_bottom/u_ice
548                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
549                  zTauB = - zTauU_ib(ji,jj) / zvel
550                  !
551                  !                 !--- Coriolis at U-points (energy conserving formulation)
552                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
553                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
554                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
555                  !
556                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
557                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
558                  !
559                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
560                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
561                  !
562                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
563                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
564                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
565                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
566                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
567                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
568                     &            ) * zmaskU(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                     &                                     + zTauE + 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_lfrelax )  &  ! static friction => slow decrease to v=0
574                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
575                     &            ) * zmaskU(ji,jj)
576                  ENDIF
577               END DO
578            END DO
579            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
580            !
581#if defined key_agrif
582!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
583            CALL agrif_interp_ice( 'U' )
584#endif
585            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
586            !
587         ELSE ! odd iterations
588            !
589            DO jj = 2, jpjm1
590               DO ji = fs_2, fs_jpim1
591                  !                 !--- tau_io/(u_oce - u_ice)
592                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
593                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
594                  !                 !--- Ocean-to-Ice stress
595                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
596                  !
597                  !                 !--- tau_bottom/u_ice
598                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
599                  zTauB = - zTauU_ib(ji,jj) / zvel
600                  !
601                  !                 !--- Coriolis at U-points (energy conserving formulation)
602                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
603                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
604                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
605                  !
606                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
607                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
608                  !
609                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
610                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
611                  !
612                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
613                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
614                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
615                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
616                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
617                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
618                     &            ) * zmaskU(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                     &                                     + zTauE + 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_lfrelax )  &  ! static friction => slow decrease to v=0
624                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
625                     &            ) * zmaskU(ji,jj)
626                  ENDIF
627               END DO
628            END DO
629            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
630            !
631#if defined key_agrif
632!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
633            CALL agrif_interp_ice( 'U' )
634#endif
635            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
636            !
637            DO jj = 2, jpjm1
638               DO ji = fs_2, fs_jpim1
639                  !                 !--- tau_io/(v_oce - v_ice)
640                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
641                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
642                  !                 !--- Ocean-to-Ice stress
643                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
644                  !
645                  !                 !--- tau_bottom/v_ice
646                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
647                  zTauB = - zTauV_ib(ji,jj) / zvel
648                  !
649                  !                 !--- Coriolis at v-points (energy conserving formulation)
650                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
651                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
652                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
653                  !
654                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
655                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
656                  !
657                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
658                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
659                  !
660                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
661                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
662                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
663                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
664                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
665                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
666                     &           ) * zmaskV(ji,jj)
667                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
668                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
669                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
670                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
671                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
672                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
673                     &            ) * zmaskV(ji,jj)
674                  ENDIF
675               END DO
676            END DO
677            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
678            !
679#if defined key_agrif
680!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
681            CALL agrif_interp_ice( 'V' )
682#endif
683            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
684            !
685         ENDIF
686
687!!$         IF(ln_ctl) THEN   ! Convergence test
688!!$            DO jj = 2 , jpjm1
689!!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
690!!$            END DO
691!!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
692!!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
693!!$         ENDIF
694         !
695         !                                                ! ==================== !
696      END DO                                              !  end loop over jter  !
697      !                                                   ! ==================== !
698      !
699      !------------------------------------------------------------------------------!
700      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
701      !------------------------------------------------------------------------------!
702      DO jj = 1, jpjm1
703         DO ji = 1, jpim1
704
705            ! shear at F points
706            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)   &
707               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
708               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
709
710         END DO
711      END DO           
712     
713      DO jj = 2, jpjm1
714         DO ji = 2, jpim1 ! no vector loop
715           
716            ! tension**2 at T points
717            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)   &
718               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
719               &   ) * r1_e1e2t(ji,jj)
720            zdt2 = zdt * zdt
721           
722            ! shear**2 at T points (doc eq. A16)
723            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
724               &   + 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)  &
725               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
726           
727            ! shear at T points
728            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
729
730            ! divergence at T points
731            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
732               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
733               &             ) * r1_e1e2t(ji,jj)
734           
735            ! delta at T points
736            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
737            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
738            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
739
740         END DO
741      END DO
742      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
743     
744      ! --- Store the stress tensor for the next time step --- !
745      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
746      pstress1_i (:,:) = zs1 (:,:)
747      pstress2_i (:,:) = zs2 (:,:)
748      pstress12_i(:,:) = zs12(:,:)
749      !
750
751      !------------------------------------------------------------------------------!
752      ! 5) diagnostics
753      !------------------------------------------------------------------------------!
754      DO jj = 1, jpj
755         DO ji = 1, jpi
756            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
757         END DO
758      END DO
759
760      ! --- divergence, shear and strength --- !
761      IF( iom_use('icediv') )   CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) )   ! divergence
762      IF( iom_use('iceshe') )   CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) )   ! shear
763      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
764
765      ! --- charge ellipse --- !
766      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
767         !
768         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
769         !         
770         DO jj = 2, jpjm1
771            DO ji = 2, jpim1
772               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
773                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
774                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
775
776               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
777
778               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
779
780!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
781!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
782!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
783!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
784               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
785               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
786               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
787            END DO
788         END DO
789         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
790         !
791         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
792         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
793         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
794         !
795         DEALLOCATE( zsig1 , zsig2 , zsig3 )
796      ENDIF
797     
798      ! --- SIMIP --- !
799      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
800         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
801         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
802         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
803
804         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
805            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
806            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
807            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
808         
809         DO jj = 2, jpjm1
810            DO ji = 2, jpim1
811               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
812               
813               ! Stress tensor invariants (normal and shear stress N/m)
814               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
815               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
816               
817               ! Stress terms of the momentum equation (N/m2)
818               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term
819               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
820               
821               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term
822               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
823               
824               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term
825               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
826               
827               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
828               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
829               
830               ! 2D ice mass, snow mass, area transport arrays (X, Y)
831               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
832               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
833               
834               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
835               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
836               
837               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
838               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
839               
840               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
841               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
842               
843            END DO
844         END DO
845         
846         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
847            &                zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
848            &                zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
849            &                zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
850                 
851         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi  , 'U', -1., zdiag_vtau_oi  , 'V', -1.,   &
852            &                zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
853            &                zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
854            &                zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
855         
856         IF( iom_use('normstr' ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
857         IF( iom_use('sheastr' ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
858         IF( iom_use('dssh_dx' ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
859         IF( iom_use('dssh_dy' ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
860         IF( iom_use('corstrx' ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
861         IF( iom_use('corstry' ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
862         IF( iom_use('intstrx' ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
863         IF( iom_use('intstry' ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
864         IF( iom_use('utau_oi' ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
865         IF( iom_use('vtau_oi' ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
866         IF( iom_use('xmtrpice') )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
867         IF( iom_use('ymtrpice') )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport
868         IF( iom_use('xmtrpsnw') )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
869         IF( iom_use('ymtrpsnw') )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
870         IF( iom_use('xatrp'   ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
871         IF( iom_use('yatrp'   ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
872
873         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
874            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
875            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
876            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
877
878      ENDIF
879      !
880   END SUBROUTINE ice_dyn_rhg_evp
881
882
883   SUBROUTINE rhg_evp_rst( cdrw, kt )
884      !!---------------------------------------------------------------------
885      !!                   ***  ROUTINE rhg_evp_rst  ***
886      !!                     
887      !! ** Purpose :   Read or write RHG file in restart file
888      !!
889      !! ** Method  :   use of IOM library
890      !!----------------------------------------------------------------------
891      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
892      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
893      !
894      INTEGER  ::   iter            ! local integer
895      INTEGER  ::   id1, id2, id3   ! local integers
896      !!----------------------------------------------------------------------
897      !
898      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
899         !                                   ! ---------------
900         IF( ln_rstart ) THEN                   !* Read the restart file
901            !
902            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
903            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
904            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
905            !
906            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
907               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
908               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
909               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
910            ELSE                                     ! start rheology from rest
911               IF(lwp) WRITE(numout,*)
912               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
913               stress1_i (:,:) = 0._wp
914               stress2_i (:,:) = 0._wp
915               stress12_i(:,:) = 0._wp
916            ENDIF
917         ELSE                                   !* Start from rest
918            IF(lwp) WRITE(numout,*)
919            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
920            stress1_i (:,:) = 0._wp
921            stress2_i (:,:) = 0._wp
922            stress12_i(:,:) = 0._wp
923         ENDIF
924         !
925      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
926         !                                   ! -------------------
927         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
928         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
929         !
930         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
931         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
932         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
933         !
934      ENDIF
935      !
936   END SUBROUTINE rhg_evp_rst
937
938#else
939   !!----------------------------------------------------------------------
940   !!   Default option         Empty module           NO SI3 sea-ice model
941   !!----------------------------------------------------------------------
942#endif
943
944   !!==============================================================================
945END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.