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

Last change on this file since 13497 was 13497, checked in by techene, 5 months ago

re-introduce comments that have been erased by loops transformation see #2525

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