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

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

changes in style - part6 - commits of the day

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