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

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

changes in style - part2 -

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