source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg_evp.F90 @ 11371

Last change on this file since 11371 was 11371, checked in by clem, 21 months ago

continue cleaning the ice outputs

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