source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90 @ 11878

Last change on this file since 11878 was 11878, checked in by stefryn, 21 months ago

add anisotropy output

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