source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90 @ 8512

Last change on this file since 8512 was 8512, checked in by clem, 3 years ago

changes in style - part5 - reaching the end

File size: 48.6 KB
Line 
1MODULE icerhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icerhg_evp  ***
4   !!   Ice rheology : sea ice rheology
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
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 cas
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 + possibility to use mEVP (Bouillon 2013)
13   !!----------------------------------------------------------------------
14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                      LIM-3 sea-ice model
17   !!----------------------------------------------------------------------
18   !!   ice_rhg_evp       : computes ice velocities
19   !!----------------------------------------------------------------------
20   USE phycst         ! Physical constant
21   USE par_oce        ! Ocean parameters
22   USE dom_oce        ! Ocean domain
23   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
24   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
25   USE ice            ! ice variables
26   USE icerdgrft      ! ice strength
27   !
28   USE lbclnk         ! Lateral Boundary Condition / MPP link
29   USE lib_mpp        ! MPP library
30   USE in_out_manager ! I/O manager
31   USE prtctl         ! Print control
32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
33   USE iom            !
34#if defined key_agrif
35   USE agrif_lim3_interp
36#endif
37   USE bdy_oce   , ONLY: ln_bdy 
38   USE bdyice 
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   ice_rhg_evp        ! routine called by icerhg.F90
44
45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
49   !! $Id: icerhg_evp.F90 8378 2017-07-26 13:55:59Z clem $
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
55      !!-------------------------------------------------------------------
56      !!                 ***  SUBROUTINE ice_rhg_evp  ***
57      !!                          EVP-C-grid
58      !!
59      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
60      !!  stress and sea-surface slope. Ice-ice interaction is described by
61      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
62      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
63      !!
64      !!  The points in the C-grid look like this, dear reader
65      !!
66      !!                              (ji,jj)
67      !!                                 |
68      !!                                 |
69      !!                      (ji-1,jj)  |  (ji,jj)
70      !!                             ---------   
71      !!                            |         |
72      !!                            | (ji,jj) |------(ji,jj)
73      !!                            |         |
74      !!                             ---------   
75      !!                     (ji-1,jj-1)     (ji,jj-1)
76      !!
77      !! ** Inputs  : - wind forcing (stress), oceanic currents
78      !!                ice total volume (vt_i) per unit area
79      !!                snow total volume (vt_s) per unit area
80      !!
81      !! ** Action  : - compute u_ice, v_ice : the components of the
82      !!                sea-ice velocity vector
83      !!              - compute delta_i, shear_i, divu_i, which are inputs
84      !!                of the ice thickness distribution
85      !!
86      !! ** Steps   : 1) Compute ice snow mass, ice strength
87      !!              2) Compute wind, oceanic stresses, mass terms and
88      !!                 coriolis terms of the momentum equation
89      !!              3) Solve the momentum equation (iterative procedure)
90      !!              4) Prevent high velocities if the ice is thin
91      !!              5) Recompute invariants of the strain rate tensor
92      !!                 which are inputs of the ITD, store stress
93      !!                 for the next time step
94      !!              6) Control prints of residual (convergence)
95      !!                 and charge ellipse.
96      !!                 The user should make sure that the parameters
97      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
98      !!                 on the charge ellipse for plastic flow
99      !!                 e.g. in the Canadian Archipelago
100      !!
101      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
102      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
103      !!              but this solution appears very unstable (see Kimmritz et al 2016)
104      !!
105      !! References : Hunke and Dukowicz, JPO97
106      !!              Bouillon et al., Ocean Modelling 2009
107      !!              Bouillon et al., Ocean Modelling 2013
108      !!-------------------------------------------------------------------
109      INTEGER, INTENT(in) ::   kt      ! time step
110      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   stress1_i, stress2_i, stress12_i
111      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   u_ice, v_ice, shear_i, divu_i, delta_i 
112      !!
113      INTEGER ::   ji, jj       ! dummy loop indices
114      INTEGER ::   jter         ! local integers
115
116      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
117      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
118      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
119      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
120      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
121      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
122      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars
123
124      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
125      REAL(wp) ::   zintb, zintn                                             ! dummy argument
126      REAL(wp) ::   zfac_x, zfac_y
127      REAL(wp) ::   zshear, zdum1, zdum2
128     
129      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
130      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
131      !
132      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
133      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
134      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
135      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
136      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
137      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
138      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
139     
140      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
141      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
142      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
143      REAL(wp), DIMENSION(jpi,jpj) ::   zpice                           ! array used for the calculation of ice surface slope:
144                                                                             !   ocean surface (ssh_m) if ice is not embedded
145                                                                             !   ice top surface if ice is embedded   
146      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
147      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
148
149      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
150      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
151      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
152
153      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
154      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
155      !! --- diags
156      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
157      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
158      !! --- SIMIP diags
159      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
160      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
161      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
162      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
163      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress 
169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
175      !!-------------------------------------------------------------------
176
177      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_rhg_evp: EVP sea-ice rheology'
178
179#if defined key_agrif 
180      CALL agrif_interp_lim3( 'U', 0, nn_nevp )   ! First interpolation of coarse values
181      CALL agrif_interp_lim3( 'V', 0, nn_nevp )
182#endif
183      !
184      !------------------------------------------------------------------------------!
185      ! 0) mask at F points for the ice
186      !------------------------------------------------------------------------------!
187      ! ocean/land mask
188      DO jj = 1, jpjm1
189         DO ji = 1, jpim1      ! NO vector opt.
190            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
191         END DO
192      END DO
193      CALL lbc_lnk( zfmask, 'F', 1._wp )
194
195      ! Lateral boundary conditions on velocity (modify zfmask)
196      zwf(:,:) = zfmask(:,:)
197      DO jj = 2, jpjm1
198         DO ji = fs_2, fs_jpim1   ! vector opt.
199            IF( zfmask(ji,jj) == 0._wp ) THEN
200               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) ) )
201            ENDIF
202         END DO
203      END DO
204      DO jj = 2, jpjm1
205         IF( zfmask(1,jj) == 0._wp ) THEN
206            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
207         ENDIF
208         IF( zfmask(jpi,jj) == 0._wp ) THEN
209            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
210         ENDIF
211      END DO
212      DO ji = 2, jpim1
213         IF( zfmask(ji,1) == 0._wp ) THEN
214            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
215         ENDIF
216         IF( zfmask(ji,jpj) == 0._wp ) THEN
217            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
218         ENDIF
219      END DO
220      CALL lbc_lnk( zfmask, 'F', 1._wp )
221
222      !------------------------------------------------------------------------------!
223      ! 1) define some variables and initialize arrays
224      !------------------------------------------------------------------------------!
225      zrhoco = rau0 * rn_cio 
226
227      ! ecc2: square of yield ellipse eccenticrity
228      ecc2    = rn_ecc * rn_ecc
229      z1_ecc2 = 1._wp / ecc2
230
231      ! Time step for subcycling
232      zdtevp   = rdt_ice / REAL( nn_nevp )
233      z1_dtevp = 1._wp / zdtevp
234
235      ! alpha parameters (Bouillon 2009)
236      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
237      zalph2 = zalph1 * z1_ecc2
238
239      ! alpha and beta parameters (Bouillon 2013)
240      !!zalph1 = 40.
241      !!zalph2 = 40.
242      !!zbeta  = 3000.
243      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001)
244
245      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
246      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
247
248      ! Initialise stress tensor
249      zs1 (:,:) = stress1_i (:,:) 
250      zs2 (:,:) = stress2_i (:,:)
251      zs12(:,:) = stress12_i(:,:)
252
253      ! Ice strength
254      CALL ice_rdgrft_strength
255
256      ! scale factors
257      DO jj = 2, jpjm1
258         DO ji = fs_2, fs_jpim1
259            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
260            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
261         END DO
262      END DO
263           
264      !
265      !------------------------------------------------------------------------------!
266      ! 2) Wind / ocean stress, mass terms, coriolis terms
267      !------------------------------------------------------------------------------!
268
269      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
270         !                                           
271         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
272         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
273         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
274         !
275         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
276         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
277         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
278         !
279         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
280         !
281      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
282         zpice(:,:) = ssh_m(:,:)
283      ENDIF
284
285      DO jj = 2, jpjm1
286         DO ji = fs_2, fs_jpim1
287
288            ! ice fraction at U-V points
289            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)
290            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)
291
292            ! Ice/snow mass at U-V points
293            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
294            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
295            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
296            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
297            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
298
299            ! Ocean currents at U-V points
300            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
301               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
302           
303            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
304               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
305
306            ! Coriolis at T points (m*f)
307            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
308
309            ! m/dt
310            zmU_t(ji,jj)    = zmassU * z1_dtevp
311            zmV_t(ji,jj)    = zmassV * z1_dtevp
312
313            ! Drag ice-atm.
314            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
315            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
316
317            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
318            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
319            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
320
321            ! masks
322            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
323            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
324
325            ! switches
326            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
327            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
328
329         END DO
330      END DO
331      CALL lbc_lnk( zmf, 'T', 1. )
332      !
333      !------------------------------------------------------------------------------!
334      ! 3) Solution of the momentum equation, iterative procedure
335      !------------------------------------------------------------------------------!
336      !
337      !                                               !----------------------!
338      DO jter = 1 , nn_nevp                           !    loop over jter    !
339         !                                            !----------------------!       
340         IF(ln_ctl) THEN   ! Convergence test
341            DO jj = 1, jpjm1
342               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
343               zv_ice(:,jj) = v_ice(:,jj)
344            END DO
345         ENDIF
346
347         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
348         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
349            DO ji = 1, jpim1
350
351               ! shear at F points
352               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)   &
353                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
354                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
355
356            END DO
357         END DO
358         CALL lbc_lnk( zds, 'F', 1. )
359
360         DO jj = 2, jpjm1
361            DO ji = 2, jpim1 ! no vector loop
362
363               ! shear**2 at T points (doc eq. A16)
364               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
365                  &   + 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)  &
366                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
367             
368               ! divergence at T points
369               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
370                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
371                  &    ) * r1_e1e2t(ji,jj)
372               zdiv2 = zdiv * zdiv
373               
374               ! tension at T points
375               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)   &
376                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
377                  &   ) * r1_e1e2t(ji,jj)
378               zdt2 = zdt * zdt
379               
380               ! delta at T points
381               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
382
383               ! P/delta at T points
384               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
385               
386               ! stress at T points
387               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
388               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
389             
390            END DO
391         END DO
392         CALL lbc_lnk( zp_delt, 'T', 1. )
393
394         DO jj = 1, jpjm1
395            DO ji = 1, jpim1
396
397               ! P/delta at F points
398               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) )
399               
400               ! stress at F points
401               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
402
403            END DO
404         END DO
405         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
406 
407
408         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
409         DO jj = 2, jpjm1
410            DO ji = fs_2, fs_jpim1               
411               !                   !--- U points
412               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
413                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
414                  &                    ) * r1_e2u(ji,jj)                                                                      &
415                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
416                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
417                  &                  ) * r1_e1e2u(ji,jj)
418                  !
419                  !                !--- V points
420               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
421                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
422                  &                    ) * r1_e1v(ji,jj)                                                                      &
423                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
424                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
425                  &                  ) * r1_e1e2v(ji,jj)
426                  !
427                  !                !--- u_ice at V point
428               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
429                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
430                  !
431                  !                !--- v_ice at U point
432               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
433                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
434               !
435            END DO
436         END DO
437         !
438         ! --- Computation of ice velocity --- !
439         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
440         !  Bouillon et al. 2009 (eq 34-35) => stable
441         IF( MOD(jter,2) == 0 ) THEN ! even iterations
442            !
443            DO jj = 2, jpjm1
444               DO ji = fs_2, fs_jpim1
445                     !                 !--- tau_io/(v_oce - v_ice)
446                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
447                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
448                     !                 !--- Ocean-to-Ice stress
449                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
450                     !
451                     !                 !--- tau_bottom/v_ice
452                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
453                  zTauB = - tau_icebfr(ji,jj) / zvel
454                     !
455                     !                 !--- Coriolis at V-points (energy conserving formulation)
456                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
457                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
458                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
459                     !
460                     !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
461                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
462                     !
463                     !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
464                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
465                     !
466                     !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
467                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
468                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
469                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
470                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
471                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
472                     &           ) * zmaskV(ji,jj)
473                     !
474                  ! Bouillon 2013
475                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
476                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
477                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
478                  !
479               END DO
480            END DO
481            CALL lbc_lnk( v_ice, 'V', -1. )
482            !
483#if defined key_agrif
484!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
485            CALL agrif_interp_lim3( 'V' )
486#endif
487            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
488            !
489            DO jj = 2, jpjm1
490               DO ji = fs_2, fs_jpim1
491                               
492                  ! tau_io/(u_oce - u_ice)
493                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
494                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
495
496                  ! Ocean-to-Ice stress
497                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
498
499                  ! tau_bottom/u_ice
500                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
501                  zTauB = - tau_icebfr(ji,jj) / zvel
502
503                  ! Coriolis at U-points (energy conserving formulation)
504                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
505                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
506                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
507                 
508                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
509                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_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 - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
513
514                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
515                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
516                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
517                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
518                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
519                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
520                     &           ) * zmaskU(ji,jj)
521                  ! Bouillon 2013
522                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
523                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
524                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
525               END DO
526            END DO
527            CALL lbc_lnk( u_ice, 'U', -1. )
528            !
529#if defined key_agrif
530!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
531            CALL agrif_interp_lim3( 'U' )
532#endif
533            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
534            !
535         ELSE ! odd iterations
536            !
537            DO jj = 2, jpjm1
538               DO ji = fs_2, fs_jpim1
539
540                  ! tau_io/(u_oce - u_ice)
541                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
542                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
543
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  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
549                  zTauB = - tau_icebfr(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 - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
561
562                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
563                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(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) + 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                  ! Bouillon 2013
570                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
571                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
572                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
573               END DO
574            END DO
575            CALL lbc_lnk( u_ice, 'U', -1. )
576            !
577#if defined key_agrif
578!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
579            CALL agrif_interp_lim3( 'U' )
580#endif
581            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
582            !
583            DO jj = 2, jpjm1
584               DO ji = fs_2, fs_jpim1
585         
586                  ! tau_io/(v_oce - v_ice)
587                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
588                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
589
590                  ! Ocean-to-Ice stress
591                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
592
593                  ! tau_bottom/v_ice
594                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
595                  ztauB = - tau_icebfr(ji,jj) / zvel
596                 
597                  ! Coriolis at V-points (energy conserving formulation)
598                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
599                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
600                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
601
602                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
603                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
604
605                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
606                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
607                 
608                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
609                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
610                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
611                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
612                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
613                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
614                     &           ) * zmaskV(ji,jj)
615                  ! Bouillon 2013
616                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
617                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
618                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
619               END DO
620            END DO
621            CALL lbc_lnk( v_ice, 'V', -1. )
622            !
623#if defined key_agrif
624!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
625            CALL agrif_interp_lim3( 'V' )
626#endif
627            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
628            !
629         ENDIF
630         
631         IF(ln_ctl) THEN   ! Convergence test
632            DO jj = 2 , jpjm1
633               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
634            END DO
635            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
636            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
637         ENDIF
638         !
639         !                                                ! ==================== !
640      END DO                                              !  end loop over jter  !
641      !                                                   ! ==================== !
642      !
643      !------------------------------------------------------------------------------!
644      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
645      !------------------------------------------------------------------------------!
646      DO jj = 1, jpjm1
647         DO ji = 1, jpim1
648
649            ! shear at F points
650            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)   &
651               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
652               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
653
654         END DO
655      END DO           
656      CALL lbc_lnk( zds, 'F', 1. )
657     
658      DO jj = 2, jpjm1
659         DO ji = 2, jpim1 ! no vector loop
660           
661            ! tension**2 at T points
662            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)   &
663               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
664               &   ) * r1_e1e2t(ji,jj)
665            zdt2 = zdt * zdt
666           
667            ! shear**2 at T points (doc eq. A16)
668            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
669               &   + 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)  &
670               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
671           
672            ! shear at T points
673            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
674
675            ! divergence at T points
676            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
677               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
678               &            ) * r1_e1e2t(ji,jj)
679           
680            ! delta at T points
681            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
682            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
683            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
684
685         END DO
686      END DO
687      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
688     
689      ! --- Store the stress tensor for the next time step --- !
690      stress1_i (:,:) = zs1 (:,:)
691      stress2_i (:,:) = zs2 (:,:)
692      stress12_i(:,:) = zs12(:,:)
693      !
694
695      !------------------------------------------------------------------------------!
696      ! 5) diagnostics
697      !------------------------------------------------------------------------------!
698      DO jj = 1, jpj
699         DO ji = 1, jpi
700            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
701         END DO
702      END DO
703
704      ! --- divergence, shear and strength --- !
705      IF( iom_use('idive')  )   CALL iom_put( "idive"  , divu_i(:,:)   * zswi(:,:) )   ! divergence
706      IF( iom_use('ishear') )   CALL iom_put( "ishear" , shear_i(:,:)  * zswi(:,:) )   ! shear
707      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
708
709      ! --- charge ellipse --- !
710      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
711         !
712         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
713         !         
714         DO jj = 2, jpjm1
715            DO ji = 2, jpim1
716               zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
717                  &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  &
718                  &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
719
720               zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
721
722               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
723
724!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
725!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
726!!               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
727!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
728               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
729               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                    ! shear stress
730               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
731            END DO
732         END DO
733         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
734         !
735         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
736         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
737         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
738         !
739         DEALLOCATE( zsig1 , zsig2 , zsig3 )
740      ENDIF
741     
742      ! --- SIMIP --- !
743      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
744         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
745         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
746         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
747
748         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
749            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
750            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
751            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
752         
753         DO jj = 2, jpjm1
754            DO ji = 2, jpim1
755               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
756               
757               ! Stress tensor invariants (normal and shear stress N/m)
758               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
759               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
760               
761               ! Stress terms of the momentum equation (N/m2)
762               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
763               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
764               
765               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
766               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
767               
768               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
769               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
770               
771               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
772               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
773               
774               ! 2D ice mass, snow mass, area transport arrays (X, Y)
775               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
776               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
777               
778               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
779               zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
780               
781               zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
782               zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
783               
784               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
785               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
786               
787            END DO
788         END DO
789         
790         CALL lbc_lnk_multi(   zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
791            &                  zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
792            &                  zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
793            &                  zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
794         
795         CALL lbc_lnk_multi(   zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1.    )
796         
797         CALL lbc_lnk_multi(   zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
798            &                  zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
799            &                  zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
800         
801         IF ( iom_use( 'normstr'  ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
802         IF ( iom_use( 'sheastr'  ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
803         IF ( iom_use( 'dssh_dx'  ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
804         IF ( iom_use( 'dssh_dy'  ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
805         IF ( iom_use( 'corstrx'  ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
806         IF ( iom_use( 'corstry'  ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
807         IF ( iom_use( 'intstrx'  ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
808         IF ( iom_use( 'intstry'  ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
809         IF ( iom_use( 'utau_oi'  ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
810         IF ( iom_use( 'vtau_oi'  ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
811         IF ( iom_use( 'xmtrpice' ) )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
812         IF ( iom_use( 'ymtrpice' ) )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport
813         IF ( iom_use( 'xmtrpsnw' ) )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
814         IF ( iom_use( 'ymtrpsnw' ) )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
815         IF ( iom_use( 'xatrp'    ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
816         IF ( iom_use( 'yatrp'    ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
817
818         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
819            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
820            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
821            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
822
823      ENDIF
824      !
825   END SUBROUTINE ice_rhg_evp
826
827#else
828   !!----------------------------------------------------------------------
829   !!   Default option         Empty module          NO LIM-3 sea-ice model
830   !!----------------------------------------------------------------------
831#endif
832
833   !!==============================================================================
834END MODULE icerhg_evp
Note: See TracBrowser for help on using the repository browser.