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

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

changes in style - part3 - move output into proper files and correct a bug in debug mode

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