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

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

source: NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rhg_evp.F90 @ 13470

Last change on this file since 13470 was 13470, checked in by smasson, 4 years ago

r4_trunk: second change of DO loops for routines to be merged, see #2523

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