source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90 @ 8554

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

changes in style - part6 - pure cosmetics

File size: 51.3 KB
Line 
1MODULE icedyn_rhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_evp  ***
4   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) 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'                                       ESIM sea-ice model
17   !!----------------------------------------------------------------------
18   !!   ice_dyn_rhg_evp       : computes ice velocities from EVP rheology
19   !!----------------------------------------------------------------------
20   USE phycst         ! Physical constant
21   USE dom_oce        ! Ocean domain
22   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
23   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
24   USE ice            ! sea-ice: ice variables
25   USE icedyn_rdgrft  ! sea-ice: ice strength
26   !
27   USE in_out_manager ! I/O manager
28   USE iom            ! I/O manager library
29   USE lib_mpp        ! MPP library
30   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
31   USE lbclnk         ! lateral boundary conditions (or mpp links)
32   USE prtctl         ! Print control
33   !
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_dyn_rhg_evp   ! called by icedyn_rhg.F90
44   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
50   !! $Id: icedyn_rhg_evp.F90 8378 2017-07-26 13:55:59Z clem $
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i )
56      !!-------------------------------------------------------------------
57      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
58      !!                          EVP-C-grid
59      !!
60      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
61      !!  stress and sea-surface slope. Ice-ice interaction is described by
62      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
63      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
64      !!
65      !!  The points in the C-grid look like this, dear reader
66      !!
67      !!                              (ji,jj)
68      !!                                 |
69      !!                                 |
70      !!                      (ji-1,jj)  |  (ji,jj)
71      !!                             ---------   
72      !!                            |         |
73      !!                            | (ji,jj) |------(ji,jj)
74      !!                            |         |
75      !!                             ---------   
76      !!                     (ji-1,jj-1)     (ji,jj-1)
77      !!
78      !! ** Inputs  : - wind forcing (stress), oceanic currents
79      !!                ice total volume (vt_i) per unit area
80      !!                snow total volume (vt_s) per unit area
81      !!
82      !! ** Action  : - compute u_ice, v_ice : the components of the
83      !!                sea-ice velocity vector
84      !!              - compute delta_i, shear_i, divu_i, which are inputs
85      !!                of the ice thickness distribution
86      !!
87      !! ** Steps   : 0) compute mask at F point
88      !!              1) Compute ice snow mass, ice strength
89      !!              2) Compute wind, oceanic stresses, mass terms and
90      !!                 coriolis terms of the momentum equation
91      !!              3) Solve the momentum equation (iterative procedure)
92      !!              4) Recompute delta, shear and divergence
93      !!                 (which are inputs of the ITD) & store stress
94      !!                 for the next time step
95      !!              5) Diagnostics including charge ellipse
96      !!
97      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
98      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
99      !!              but this solution appears very unstable (see Kimmritz et al 2016)
100      !!
101      !! References : Hunke and Dukowicz, JPO97
102      !!              Bouillon et al., Ocean Modelling 2009
103      !!              Bouillon et al., Ocean Modelling 2013
104      !!-------------------------------------------------------------------
105      INTEGER, INTENT(in) ::   kt      ! time step
106      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i
107      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i 
108      !!
109      INTEGER ::   ji, jj       ! dummy loop indices
110      INTEGER ::   jter         ! local integers
111
112      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
113      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
114      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
115      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
116      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
117      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
118      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars
119
120      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
121      REAL(wp) ::   zintb, zintn                                             ! dummy argument
122      REAL(wp) ::   zfac_x, zfac_y
123      REAL(wp) ::   zshear, zdum1, zdum2
124     
125      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
126      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
127      !
128      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
129      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
130      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
131      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
132      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
133      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
134      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
135     
136      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
137      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
138      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
139      REAL(wp), DIMENSION(jpi,jpj) ::   zpice                           ! array used for the calculation of ice surface slope:
140                                                                             !   ocean surface (ssh_m) if ice is not embedded
141                                                                             !   ice top surface if ice is embedded   
142      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
143      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
144
145      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
146      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
147      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
148
149      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
150      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
151      !! --- diags
152      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
153      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
154      !! --- SIMIP diags
155      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
156      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
157      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
158      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
159      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
160      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
161      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
162      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
163      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress 
165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
171      !!-------------------------------------------------------------------
172
173      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
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 (:,:) = pstress1_i (:,:) 
246      zs2 (:,:) = pstress2_i (:,:)
247      zs12(:,:) = pstress12_i(:,:)
248
249      ! Ice strength
250      CALL ice_strength
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) = pu_ice(:,jj) ! velocity at previous time step
339               zv_ice(:,jj) = pv_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) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
349                  &         + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_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) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj)   &
366                  &    + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1)   &
367                  &    ) * r1_e1e2t(ji,jj)
368               zdiv2 = zdiv * zdiv
369               
370               ! tension at T points
371               zdt  = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
372                  &   - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_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 * ( ( pu_ice(ji,jj  ) + pu_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
425                  &                     + ( pu_ice(ji,jj+1) + pu_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 * ( ( pv_ice(ji  ,jj) + pv_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
429                  &                     + ( pv_ice(ji+1,jj) + pv_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( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_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) - pv_ice(ji,jj) )
446                     !
447                     !                 !--- tau_bottom/v_ice
448                  zvel  = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_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  ) * pu_ice(ji,jj  ) + e2u(ji-1,jj  ) * pu_ice(ji-1,jj  ) )  &
454                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_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                  pv_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * pv_ice(ji,jj)                            &  ! previous velocity
464                     &                                      + zTauE + zTauO * pv_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 ) * pv_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                  !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_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( pv_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( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_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) - pu_ice(ji,jj) )
494
495                  ! tau_bottom/u_ice
496                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_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) * pv_ice(ji  ,jj) + e1v(ji  ,jj-1) * pv_ice(ji  ,jj-1) )  &
502                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_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                  pu_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * pu_ice(ji,jj)                            &  ! previous velocity
512                     &                                      + zTauE + zTauO * pu_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 ) * pu_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                  !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_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( pu_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( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_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) - pu_ice(ji,jj) )
542
543                  ! tau_bottom/u_ice
544                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_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) * pv_ice(ji  ,jj) + e1v(ji  ,jj-1) * pv_ice(ji  ,jj-1) )  &
550                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_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                  pu_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * pu_ice(ji,jj)                            &  ! previous velocity
560                     &                                      + zTauE + zTauO * pu_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 ) * pu_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                  !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_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( pu_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( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_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) - pv_ice(ji,jj) )
588
589                  ! tau_bottom/v_ice
590                  zvel  = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_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  ) * pu_ice(ji,jj  ) + e2u(ji-1,jj  ) * pu_ice(ji-1,jj  ) )  &
596                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_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                  pv_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * pv_ice(ji,jj)                            &  ! previous velocity
606                     &                                      + zTauE + zTauO * pv_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 ) * pv_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                  !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_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( pv_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( pu_ice(:,jj) - zu_ice(:,jj) ), ABS( pv_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) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
647               &         + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_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  = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
659               &   - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_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            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
670
671            ! divergence at T points
672            pdivu_i(ji,jj) = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj)   &
673               &             + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1)   &
674               &             ) * r1_e1e2t(ji,jj)
675           
676            ! delta at T points
677            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
678            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
679            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
680
681         END DO
682      END DO
683      CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
684     
685      ! --- Store the stress tensor for the next time step --- !
686      pstress1_i (:,:) = zs1 (:,:)
687      pstress2_i (:,:) = zs2 (:,:)
688      pstress12_i(:,:) = zs12(:,:)
689      !
690
691      !------------------------------------------------------------------------------!
692      ! 5) diagnostics
693      !------------------------------------------------------------------------------!
694      DO jj = 1, jpj
695         DO ji = 1, jpi
696            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
697         END DO
698      END DO
699
700      ! --- divergence, shear and strength --- !
701      IF( iom_use('idive')  )   CALL iom_put( "idive"  , pdivu_i(:,:)   * zswi(:,:) )   ! divergence
702      IF( iom_use('ishear') )   CALL iom_put( "ishear" , pshear_i(:,:)  * zswi(:,:) )   ! shear
703      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
704
705      ! --- charge ellipse --- !
706      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
707         !
708         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
709         !         
710         DO jj = 2, jpjm1
711            DO ji = 2, jpim1
712               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
713                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
714                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
715
716               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
717
718               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
719
720!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
721!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
722!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
723!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
724               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
725               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                    ! shear stress
726               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
727            END DO
728         END DO
729         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
730         !
731         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
732         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
733         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
734         !
735         DEALLOCATE( zsig1 , zsig2 , zsig3 )
736      ENDIF
737     
738      ! --- SIMIP --- !
739      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
740         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
741         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
742         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
743
744         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
745            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
746            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
747            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
748         
749         DO jj = 2, jpjm1
750            DO ji = 2, jpim1
751               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
752               
753               ! Stress tensor invariants (normal and shear stress N/m)
754               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
755               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
756               
757               ! Stress terms of the momentum equation (N/m2)
758               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
759               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
760               
761               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
762               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
763               
764               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
765               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
766               
767               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
768               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
769               
770               ! 2D ice mass, snow mass, area transport arrays (X, Y)
771               zfac_x = 0.5 * pu_ice(ji,jj) * e2u(ji,jj) * rswitch
772               zfac_y = 0.5 * pv_ice(ji,jj) * e1v(ji,jj) * rswitch
773               
774               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
775               zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
776               
777               zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
778               zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
779               
780               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
781               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
782               
783            END DO
784         END DO
785         
786         CALL lbc_lnk_multi(   zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
787            &                  zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
788            &                  zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
789            &                  zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
790         
791         CALL lbc_lnk_multi(   zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1.    )
792         
793         CALL lbc_lnk_multi(   zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
794            &                  zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
795            &                  zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
796         
797         IF ( iom_use( 'normstr'  ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
798         IF ( iom_use( 'sheastr'  ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
799         IF ( iom_use( 'dssh_dx'  ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
800         IF ( iom_use( 'dssh_dy'  ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
801         IF ( iom_use( 'corstrx'  ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
802         IF ( iom_use( 'corstry'  ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
803         IF ( iom_use( 'intstrx'  ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
804         IF ( iom_use( 'intstry'  ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
805         IF ( iom_use( 'utau_oi'  ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
806         IF ( iom_use( 'vtau_oi'  ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
807         IF ( iom_use( 'xmtrpice' ) )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
808         IF ( iom_use( 'ymtrpice' ) )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport
809         IF ( iom_use( 'xmtrpsnw' ) )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
810         IF ( iom_use( 'ymtrpsnw' ) )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
811         IF ( iom_use( 'xatrp'    ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
812         IF ( iom_use( 'yatrp'    ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
813
814         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
815            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
816            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
817            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
818
819      ENDIF
820      !
821   END SUBROUTINE ice_dyn_rhg_evp
822
823   SUBROUTINE rhg_evp_rst( cdrw, kt )
824      !!---------------------------------------------------------------------
825      !!                   ***  ROUTINE rhg_evp_rst  ***
826      !!                     
827      !! ** Purpose :   Read or write RHG file in restart file
828      !!
829      !! ** Method  :   use of IOM library
830      !!----------------------------------------------------------------------
831      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
832      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
833      !
834      INTEGER  ::   iter            ! local integer
835      INTEGER  ::   id1, id2, id3   ! local integers
836      !!----------------------------------------------------------------------
837      !
838      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
839         !                                   ! ---------------
840         IF( ln_rstart ) THEN                   !* Read the restart file
841            !
842            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
843            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
844            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
845            !
846            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
847               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
848               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
849               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
850            ELSE                                     ! start rheology from rest
851               IF(lwp) WRITE(numout,*) '   ==>>   previous run without rheology, set stresses to 0'
852               stress1_i (:,:) = 0._wp
853               stress2_i (:,:) = 0._wp
854               stress12_i(:,:) = 0._wp
855            ENDIF
856         ELSE                                   !* Start from rest
857            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set stresses to 0'
858            stress1_i (:,:) = 0._wp
859            stress2_i (:,:) = 0._wp
860            stress12_i(:,:) = 0._wp
861         ENDIF
862         !
863      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
864         !                                   ! -------------------
865         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
866         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
867         !
868         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
869         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
870         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
871         !
872      ENDIF
873      !
874   END SUBROUTINE rhg_evp_rst
875
876#else
877   !!----------------------------------------------------------------------
878   !!   Default option         Empty module           NO ESIM sea-ice model
879   !!----------------------------------------------------------------------
880#endif
881
882   !!==============================================================================
883END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.