source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg_evp.F90 @ 10365

Last change on this file since 10365 was 10365, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with dev_r9866_HPC_03_globcom, see #2133

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