source: NEMO/branches/UKMO/dev_r10037_GPU/src/ICE/icedyn_rhg_evp.F90 @ 11888

Last change on this file since 11888 was 11888, checked in by andmirek, 19 months ago

Ticket #2197, small changes to improve code structure

File size: 56.5 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) adaptation to new model
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 case
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 + mEVP (Bouillon 2013)
13   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
14   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
15   !!----------------------------------------------------------------------
16#if defined key_si3
17   !!----------------------------------------------------------------------
18   !!   'key_si3'                                       SI3 sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
21   !!   rhg_evp_rst     : read/write EVP fields in ice restart
22   !!----------------------------------------------------------------------
23   USE phycst         ! Physical constant
24   USE dom_oce        ! Ocean domain
25   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
26   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
27   USE ice            ! sea-ice: ice variables
28   USE icedyn_rdgrft  ! sea-ice: ice strength
29   USE bdy_oce , ONLY : ln_bdy 
30   USE bdyice 
31#if defined key_agrif
32   USE agrif_ice_interp
33#endif
34   !
35   USE in_out_manager ! I/O manager
36   USE iom            ! I/O manager library
37   USE lib_mpp        ! MPP library
38   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
39   USE lbclnk         ! lateral boundary conditions (or mpp links)
40   USE prtctl         ! Print control
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
46   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
47
48   !! * Substitutions
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
58      !!-------------------------------------------------------------------
59      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
60      !!                             EVP-C-grid
61      !!
62      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
63      !!  stress and sea-surface slope. Ice-ice interaction is described by
64      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
65      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
66      !!
67      !!  The points in the C-grid look like this, dear reader
68      !!
69      !!                              (ji,jj)
70      !!                                 |
71      !!                                 |
72      !!                      (ji-1,jj)  |  (ji,jj)
73      !!                             ---------   
74      !!                            |         |
75      !!                            | (ji,jj) |------(ji,jj)
76      !!                            |         |
77      !!                             ---------   
78      !!                     (ji-1,jj-1)     (ji,jj-1)
79      !!
80      !! ** Inputs  : - wind forcing (stress), oceanic currents
81      !!                ice total volume (vt_i) per unit area
82      !!                snow total volume (vt_s) per unit area
83      !!
84      !! ** Action  : - compute u_ice, v_ice : the components of the
85      !!                sea-ice velocity vector
86      !!              - compute delta_i, shear_i, divu_i, which are inputs
87      !!                of the ice thickness distribution
88      !!
89      !! ** Steps   : 0) compute mask at F point
90      !!              1) Compute ice snow mass, ice strength
91      !!              2) Compute wind, oceanic stresses, mass terms and
92      !!                 coriolis terms of the momentum equation
93      !!              3) Solve the momentum equation (iterative procedure)
94      !!              4) Recompute delta, shear and divergence
95      !!                 (which are inputs of the ITD) & store stress
96      !!                 for the next time step
97      !!              5) Diagnostics including charge ellipse
98      !!
99      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
100      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
101      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
102      !!              (i.e. more stable and better convergence)
103      !!
104      !! References : Hunke and Dukowicz, JPO97
105      !!              Bouillon et al., Ocean Modelling 2009
106      !!              Bouillon et al., Ocean Modelling 2013
107      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
108      !!-------------------------------------------------------------------
109      USE scoce, ONLY : z1_e1t0 => scr2D1, z1_e2t0 => scr2D2, &                                       ! scale factors
110                        zp_delt => scr2D3, &                                                          ! P/delta at T points
111                        zbeta => scr2D4, &                                                            ! beta coef from Kimmritz 2017
112      !
113                        zdt_m => scr2D5, &                                                            ! (dt / ice-snow_mass) on T points
114                        zaU => scr2D6  , zaV => scr2D7, &                                             ! ice fraction on U/V points
115                        zmU_t => scr2D8, zmV_t => scr2D9, &                                           ! (ice-snow_mass / dt) on U/V points
116                        zmf => scr2D10, &                                                             ! coriolis parameter at T points
117                        zTauU_ia => scr2D11, ztauV_ia => scr2D12, &                                   ! ice-atm. stress at U-V points
118                        zspgU => scr2D13, zspgV => scr2D14, &                                         ! surface pressure gradient at U/V points
119                        v_oceU => scr2D15, u_oceV => scr2D16, v_iceU => scr2D17, u_iceV => scr2D18, & ! ocean/ice u/v component on V/U points                           
120                        zfU => scr2D19  , zfV => scr2D20, &                                           ! internal stresses
121      !
122                        zds => scr2D21, &                                                             ! shear
123                        zs1 => scr2D22, zs2 => scr2D23, zs12 => scr2D38, &                            ! stress tensor components
124                        zu_ice => scr2D24, zv_ice => scr2D25, zresr => scr2D26, &                     ! check convergence
125                        zpice  => scr2D27, &                                                          ! array used for the calculation of ice surface slope:
126      !                                                                                               !    ocean surface (ssh_m) if ice is not embedded
127      !                                                                                               !    ice top surface if ice is embedded   
128                        zCorx => scr2D28, zCory => scr2D29, &                                         ! Coriolis stress array
129                        ztaux_oi => scr2D30, ztauy_oi => scr2D31, &                                   ! Ocean-to-ice stress array
130      !
131                        zswitchU => scr2D32, zswitchV => scr2D33, &                                   ! dummy arrays
132                        zmaskU => scr2D34, zmaskV => scr2D35, &                                       ! mask for ice presence
133                        zfmask => scr2D36, zwf => scr2D37                                             ! mask at F points for the ice
134
135                       
136
137      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
138      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
139      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
140      !!
141      INTEGER ::   ji, jj       ! dummy loop indices
142      INTEGER ::   jter         ! local integers
143      !
144      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
145      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
146      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
147      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
148      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass
149      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
150      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars
151      !
152      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
153      REAL(wp) ::   zintb, zintn                                        ! dummy argument
154      REAL(wp) ::   zfac_x, zfac_y
155      REAL(wp) ::   zshear, zdum1, zdum2
156      !
157      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
158      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
159      !! --- diags
160      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
161      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
162      !! --- SIMIP diags
163      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress 
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
179      !!-------------------------------------------------------------------
180
181      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
182      !
183!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
184      !------------------------------------------------------------------------------!
185      ! 0) mask at F points for the ice
186      !------------------------------------------------------------------------------!
187      ! ocean/land mask
188      DO jj = 1, jpjm1
189         DO ji = 1, jpim1      ! NO vector opt.
190            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
191         END DO
192      END DO
193      CALL lbc_lnk( zfmask, 'F', 1._wp )
194
195      ! Lateral boundary conditions on velocity (modify zfmask)
196      zwf(:,:) = zfmask(:,:)
197      DO jj = 2, jpjm1
198         DO ji = fs_2, fs_jpim1   ! vector opt.
199            IF( zfmask(ji,jj) == 0._wp ) THEN
200               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
201            ENDIF
202         END DO
203      END DO
204      DO jj = 2, jpjm1
205         IF( zfmask(1,jj) == 0._wp ) THEN
206            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
207         ENDIF
208         IF( zfmask(jpi,jj) == 0._wp ) THEN
209            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
210         ENDIF
211      END DO
212      DO ji = 2, jpim1
213         IF( zfmask(ji,1) == 0._wp ) THEN
214            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
215         ENDIF
216         IF( zfmask(ji,jpj) == 0._wp ) THEN
217            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
218         ENDIF
219      END DO
220      CALL lbc_lnk( zfmask, 'F', 1._wp )
221
222      !------------------------------------------------------------------------------!
223      ! 1) define some variables and initialize arrays
224      !------------------------------------------------------------------------------!
225      zrhoco = rau0 * rn_cio 
226
227      ! ecc2: square of yield ellipse eccenticrity
228      ecc2    = rn_ecc * rn_ecc
229      z1_ecc2 = 1._wp / ecc2
230
231      ! Time step for subcycling
232      zdtevp   = rdt_ice / REAL( nn_nevp )
233      z1_dtevp = 1._wp / zdtevp
234
235      ! alpha parameters (Bouillon 2009)
236      IF( .NOT. ln_aEVP ) THEN
237         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
238         zalph2 = zalph1 * z1_ecc2
239
240         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
241         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
242      ENDIF
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 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
290            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
291            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * 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            ! dt/m at T points (for alpha and beta coefficients)
306            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
307           
308            ! m/dt
309            zmU_t(ji,jj)    = zmassU * z1_dtevp
310            zmV_t(ji,jj)    = zmassV * z1_dtevp
311           
312            ! Drag ice-atm.
313            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
314            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
315
316            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
317            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
318            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
319
320            ! masks
321            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
322            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
323
324            ! switches
325            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
326            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
327
328         END DO
329      END DO
330      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )
331      !
332      !------------------------------------------------------------------------------!
333      ! 3) Solution of the momentum equation, iterative procedure
334      !------------------------------------------------------------------------------!
335      !
336      !                                               !----------------------!
337      DO jter = 1 , nn_nevp                           !    loop over jter    !
338         !                                            !----------------------!       
339         IF(ln_ctl) THEN   ! Convergence test
340            DO jj = 1, jpjm1
341               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
342               zv_ice(:,jj) = v_ice(:,jj)
343            END DO
344         ENDIF
345
346         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
347         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
348            DO ji = 1, jpim1
349
350               ! shear at F points
351               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
352                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
353                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
354
355            END DO
356         END DO
357         CALL lbc_lnk( zds, 'F', 1. )
358
359         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
360            DO ji = 2, jpi ! no vector loop
361
362               ! shear**2 at T points (doc eq. A16)
363               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
364                  &   + 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)  &
365                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
366             
367               ! divergence at T points
368               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
369                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
370                  &    ) * r1_e1e2t(ji,jj)
371               zdiv2 = zdiv * zdiv
372               
373               ! tension at T points
374               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
375                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
376                  &   ) * r1_e1e2t(ji,jj)
377               zdt2 = zdt * zdt
378               
379               ! delta at T points
380               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
381
382               ! P/delta at T points
383               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
384
385               ! alpha & beta for aEVP
386               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
387               !   alpha = beta = sqrt(4*gamma)
388               IF( ln_aEVP ) THEN
389                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
390                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
391                  zalph2   = zalph1
392                  z1_alph2 = z1_alph1
393               ENDIF
394               
395               ! stress at T points
396               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
397               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
398             
399            END DO
400         END DO
401         CALL lbc_lnk( zp_delt, 'T', 1. )
402
403         IF( ln_aEVP ) THEN
404            DO jj = 1, jpjm1
405               DO ji = 1, jpim1
406
407                  ! alpha & beta for aEVP
408                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
409                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
410                  zbeta(ji,jj) = zalph2
411               
412                  ! P/delta at F points
413                  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) )
414               
415                  ! stress at F points
416                  zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
417
418               END DO
419            END DO
420         ELSE
421            DO jj = 1, jpjm1
422               DO ji = 1, jpim1
423
424               ! P/delta at F points
425               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) )
426
427               ! stress at F points
428               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
429
430               END DO
431            END DO
432         ENDIF
433
434         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
435         DO jj = 2, jpjm1
436            DO ji = fs_2, fs_jpim1               
437               !                   !--- U points
438               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
439                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
440                  &                    ) * r1_e2u(ji,jj)                                                                      &
441                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
442                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
443                  &                  ) * r1_e1e2u(ji,jj)
444               !
445               !                !--- V points
446               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
447                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
448                  &                    ) * r1_e1v(ji,jj)                                                                      &
449                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
450                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
451                  &                  ) * r1_e1e2v(ji,jj)
452               !
453               !                !--- u_ice at V point
454               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
455                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
456               !
457               !                !--- v_ice at U point
458               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
459                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
460            END DO
461         END DO
462         !
463         ! --- Computation of ice velocity --- !
464         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
465         !  Bouillon et al. 2009 (eq 34-35) => stable
466         IF( MOD(jter,2) == 0 ) THEN ! even iterations
467            !
468            DO jj = 2, jpjm1
469               DO ji = fs_2, fs_jpim1
470                  !                 !--- tau_io/(v_oce - v_ice)
471                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
472                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
473                  !                 !--- Ocean-to-Ice stress
474                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
475                  !
476                  !                 !--- tau_bottom/v_ice
477                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
478                  zTauB = - tau_icebfr(ji,jj) / zvel
479                  !
480                  !                 !--- Coriolis at V-points (energy conserving formulation)
481                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
482                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
483                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
484                  !
485                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
486                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
487                  !
488                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
489                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
490                  !
491                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
492                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
493                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
494                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
495                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
496                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
497                     &           ) * zmaskV(ji,jj)
498                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
499                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
500                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
501                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
502                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
503                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
504                     &            ) * zmaskV(ji,jj)
505                  ENDIF
506               END DO
507            END DO
508            CALL lbc_lnk( v_ice, 'V', -1. )
509            !
510#if defined key_agrif
511!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
512            CALL agrif_interp_ice( 'V' )
513#endif
514            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
515            !
516            DO jj = 2, jpjm1
517               DO ji = fs_2, fs_jpim1         
518                  !                 !--- tau_io/(u_oce - u_ice)
519                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
520                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
521                  !                 !--- Ocean-to-Ice stress
522                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
523                  !
524                  !                 !--- tau_bottom/u_ice
525                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
526                  zTauB = - tau_icebfr(ji,jj) / zvel
527                  !
528                  !                 !--- Coriolis at U-points (energy conserving formulation)
529                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
530                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
531                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
532                  !
533                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
534                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
535                  !
536                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
537                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
538                  !
539                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
540                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
541                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
542                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
543                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
544                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
545                     &            ) * zmaskU(ji,jj)
546                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
547                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
548                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
549                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
550                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
551                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
552                     &            ) * zmaskU(ji,jj)
553                  ENDIF
554               END DO
555            END DO
556            CALL lbc_lnk( u_ice, 'U', -1. )
557            !
558#if defined key_agrif
559!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
560            CALL agrif_interp_ice( 'U' )
561#endif
562            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
563            !
564         ELSE ! odd iterations
565            !
566            DO jj = 2, jpjm1
567               DO ji = fs_2, fs_jpim1
568                  !                 !--- tau_io/(u_oce - u_ice)
569                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
570                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
571                  !                 !--- Ocean-to-Ice stress
572                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
573                  !
574                  !                 !--- tau_bottom/u_ice
575                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
576                  zTauB = - tau_icebfr(ji,jj) / zvel
577                  !
578                  !                 !--- Coriolis at U-points (energy conserving formulation)
579                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
580                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
581                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
582                  !
583                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
584                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
585                  !
586                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
587                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
588                  !
589                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
590                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
591                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
592                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
593                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
594                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
595                     &            ) * zmaskU(ji,jj)
596                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
597                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
598                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
599                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
600                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
601                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
602                     &            ) * zmaskU(ji,jj)
603                  ENDIF
604               END DO
605            END DO
606            CALL lbc_lnk( u_ice, 'U', -1. )
607            !
608#if defined key_agrif
609!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
610            CALL agrif_interp_ice( 'U' )
611#endif
612            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
613            !
614            DO jj = 2, jpjm1
615               DO ji = fs_2, fs_jpim1
616                  !                 !--- tau_io/(v_oce - v_ice)
617                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
618                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
619                  !                 !--- Ocean-to-Ice stress
620                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
621                  !
622                  !                 !--- tau_bottom/v_ice
623                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
624                  ztauB = - tau_icebfr(ji,jj) / zvel
625                  !
626                  !                 !--- Coriolis at v-points (energy conserving formulation)
627                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
628                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
629                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
630                  !
631                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
632                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
633                  !
634                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
635                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
636                  !
637                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
638                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
639                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
640                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
641                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
642                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
643                     &           ) * zmaskV(ji,jj)
644                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
645                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
646                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
647                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
648                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
649                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
650                     &            ) * zmaskV(ji,jj)
651                  ENDIF
652               END DO
653            END DO
654            CALL lbc_lnk( v_ice, 'V', -1. )
655            !
656#if defined key_agrif
657!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
658            CALL agrif_interp_ice( 'V' )
659#endif
660            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
661            !
662         ENDIF
663         
664         IF(ln_ctl) THEN   ! Convergence test
665            DO jj = 2 , jpjm1
666               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
667            END DO
668            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
669            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
670         ENDIF
671         !
672         !                                                ! ==================== !
673      END DO                                              !  end loop over jter  !
674      !                                                   ! ==================== !
675      !
676      !------------------------------------------------------------------------------!
677      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
678      !------------------------------------------------------------------------------!
679      DO jj = 1, jpjm1
680         DO ji = 1, jpim1
681
682            ! shear at F points
683            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
684               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
685               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
686
687         END DO
688      END DO           
689     
690      DO jj = 2, jpjm1
691         DO ji = 2, jpim1 ! no vector loop
692           
693            ! tension**2 at T points
694            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
695               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
696               &   ) * r1_e1e2t(ji,jj)
697            zdt2 = zdt * zdt
698           
699            ! shear**2 at T points (doc eq. A16)
700            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
701               &   + 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)  &
702               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
703           
704            ! shear at T points
705            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
706
707            ! divergence at T points
708            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
709               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
710               &             ) * r1_e1e2t(ji,jj)
711           
712            ! delta at T points
713            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
714            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
715            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
716
717         END DO
718      END DO
719      CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
720     
721      ! --- Store the stress tensor for the next time step --- !
722      CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
723      pstress1_i (:,:) = zs1 (:,:)
724      pstress2_i (:,:) = zs2 (:,:)
725      pstress12_i(:,:) = zs12(:,:)
726      !
727
728      !------------------------------------------------------------------------------!
729      ! 5) diagnostics
730      !------------------------------------------------------------------------------!
731      DO jj = 1, jpj
732         DO ji = 1, jpi
733            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
734         END DO
735      END DO
736
737      ! --- divergence, shear and strength --- !
738      IF( iom_use('icediv') )   CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) )   ! divergence
739      IF( iom_use('iceshe') )   CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) )   ! shear
740      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
741
742      ! --- charge ellipse --- !
743      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
744         !
745         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
746         !         
747         DO jj = 2, jpjm1
748            DO ji = 2, jpim1
749               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
750                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
751                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
752
753               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
754
755               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
756
757!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
758!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
759!!               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
760!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
761               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
762               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
763               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
764            END DO
765         END DO
766         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
767         !
768         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
769         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
770         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
771         !
772         DEALLOCATE( zsig1 , zsig2 , zsig3 )
773      ENDIF
774     
775      ! --- SIMIP --- !
776      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
777         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
778         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
779         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
780
781         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
782            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
783            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
784            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
785         
786         DO jj = 2, jpjm1
787            DO ji = 2, jpim1
788               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
789               
790               ! Stress tensor invariants (normal and shear stress N/m)
791               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
792               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
793               
794               ! Stress terms of the momentum equation (N/m2)
795               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term
796               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
797               
798               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term
799               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
800               
801               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term
802               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
803               
804               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
805               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
806               
807               ! 2D ice mass, snow mass, area transport arrays (X, Y)
808               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
809               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
810               
811               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
812               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
813               
814               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
815               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
816               
817               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
818               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
819               
820            END DO
821         END DO
822         
823         CALL lbc_lnk_multi( zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
824            &                zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
825            &                zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
826            &                zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
827                 
828         CALL lbc_lnk_multi( zdiag_utau_oi  , 'U', -1., zdiag_vtau_oi  , 'V', -1.,   &
829            &                zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
830            &                zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
831            &                zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
832         
833         IF( iom_use('normstr' ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
834         IF( iom_use('sheastr' ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
835         IF( iom_use('dssh_dx' ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
836         IF( iom_use('dssh_dy' ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
837         IF( iom_use('corstrx' ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
838         IF( iom_use('corstry' ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
839         IF( iom_use('intstrx' ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
840         IF( iom_use('intstry' ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
841         IF( iom_use('utau_oi' ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
842         IF( iom_use('vtau_oi' ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
843         IF( iom_use('xmtrpice') )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
844         IF( iom_use('ymtrpice') )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport
845         IF( iom_use('xmtrpsnw') )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
846         IF( iom_use('ymtrpsnw') )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
847         IF( iom_use('xatrp'   ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
848         IF( iom_use('yatrp'   ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
849
850         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
851            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
852            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
853            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
854
855      ENDIF
856      !
857   END SUBROUTINE ice_dyn_rhg_evp
858
859
860   SUBROUTINE rhg_evp_rst( cdrw, kt )
861      !!---------------------------------------------------------------------
862      !!                   ***  ROUTINE rhg_evp_rst  ***
863      !!                     
864      !! ** Purpose :   Read or write RHG file in restart file
865      !!
866      !! ** Method  :   use of IOM library
867      !!----------------------------------------------------------------------
868      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
869      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
870      !
871      INTEGER  ::   iter            ! local integer
872      INTEGER  ::   id1, id2, id3   ! local integers
873      !!----------------------------------------------------------------------
874      !
875      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
876         !                                   ! ---------------
877         IF( ln_rstart ) THEN                   !* Read the restart file
878            !
879            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
880            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
881            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
882            !
883            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
884               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
885               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
886               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
887            ELSE                                     ! start rheology from rest
888               IF(lwp) WRITE(numout,*)
889               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
890               stress1_i (:,:) = 0._wp
891               stress2_i (:,:) = 0._wp
892               stress12_i(:,:) = 0._wp
893            ENDIF
894         ELSE                                   !* Start from rest
895            IF(lwp) WRITE(numout,*)
896            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
897            stress1_i (:,:) = 0._wp
898            stress2_i (:,:) = 0._wp
899            stress12_i(:,:) = 0._wp
900         ENDIF
901         !
902      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
903         !                                   ! -------------------
904         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
905         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
906         !
907         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
908         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
909         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
910         !
911      ENDIF
912      !
913   END SUBROUTINE rhg_evp_rst
914
915#else
916   !!----------------------------------------------------------------------
917   !!   Default option         Empty module           NO SI3 sea-ice model
918   !!----------------------------------------------------------------------
919#endif
920
921   !!==============================================================================
922END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.