New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rhg_evp.F90 in NEMO/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 @ 9604

Last change on this file since 9604 was 9604, checked in by clem, 7 years ago

change history of the ice routines

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