source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icedyn_rhg_evp.F90 @ 12475

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