ecearth3/sources/nemo-3.6/NEMO/LIM_SRC_3/limrhg.F90

MODULE limrhg
   !!======================================================================
   !!                     ***  MODULE  limrhg  ***
   !!   Ice rheology : sea ice rheology
   !!======================================================================
   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas
   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
   !!            3.5  !  2012-08  (R. Benshila)  AGRIF 
   !!            3.6  !  2016-06  (C. Rousset) Rewriting (conserves energy)
   !!----------------------------------------------------------------------
#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
   !!----------------------------------------------------------------------
   !!   'key_lim3'               OR                     LIM-3 sea-ice model
   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model
   !!----------------------------------------------------------------------
   !!   lim_rhg       : computes ice velocities
   !!----------------------------------------------------------------------
   USE phycst         ! Physical constant
   USE oce     , ONLY :  snwice_mass, snwice_mass_b
   USE par_oce        ! Ocean parameters
   USE dom_oce        ! Ocean domain
   USE sbc_oce        ! Surface boundary condition: ocean fields
   USE sbc_ice        ! Surface boundary condition: ice fields
#if defined key_lim3
   USE ice            ! LIM-3: ice variables
   USE dom_ice        ! LIM-3: ice domain
   USE limitd_me      ! LIM-3: 
#else
   USE ice_2          ! LIM-2: ice variables
   USE dom_ice_2      ! LIM-2: ice domain
#endif
   USE lbclnk         ! Lateral Boundary Condition / MPP link
   USE lib_mpp        ! MPP library
   USE wrk_nemo       ! work arrays
   USE in_out_manager ! I/O manager
   USE prtctl         ! Print control
   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
#if defined key_agrif && defined key_lim2
   USE agrif_lim2_interp
#endif
#if defined key_bdy
   USE bdyice_lim
#endif

   IMPLICIT NONE
   PRIVATE

   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)

   !! * Substitutions
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
   !! $Id: limrhg.F90 4990 2014-12-15 16:42:49Z timgraham $
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE lim_rhg( k_j1, k_jpj )
      !!-------------------------------------------------------------------
      !!                 ***  SUBROUTINE lim_rhg  ***
      !!                          EVP-C-grid
      !!
      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
      !!  stress and sea-surface slope. Ice-ice interaction is described by 
      !!  a non-linear elasto-viscous-plastic (EVP) law including shear 
      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).	
      !!
      !!  The points in the C-grid look like this, dear reader
      !!
      !!                              (ji,jj)
      !!                                 |
      !!                                 |
      !!                      (ji-1,jj)  |  (ji,jj)
      !!                             ---------   
      !!                            |         |
      !!                            | (ji,jj) |------(ji,jj)
      !!                            |         |
      !!                             ---------   
      !!                     (ji-1,jj-1)     (ji,jj-1)
      !!
      !! ** Inputs  : - wind forcing (stress), oceanic currents
      !!                ice total volume (vt_i) per unit area
      !!                snow total volume (vt_s) per unit area
      !!
      !! ** Action  : - compute u_ice, v_ice : the components of the 
      !!                sea-ice velocity vector
      !!              - compute delta_i, shear_i, divu_i, which are inputs
      !!                of the ice thickness distribution
      !!
      !! ** Steps   : 1) Compute ice snow mass, ice strength 
      !!              2) Compute wind, oceanic stresses, mass terms and
      !!                 coriolis terms of the momentum equation
      !!              3) Solve the momentum equation (iterative procedure)
      !!              4) Recompute invariants of the strain rate tensor
      !!                 which are inputs of the ITD, store stress
      !!                 for the next time step
      !!              5) Control prints of residual (convergence)
      !!                 and charge ellipse.
      !!                 The user should make sure that the parameters
      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
      !!                 on the charge ellipse for plastic flow
      !!                 e.g. in the Canadian Archipelago
      !!
      !! ** Notes   : Boundary condition for ice is chosen no-slip 
      !!              but can be adjusted with param rn_shlat
      !!
      !! References : Hunke and Dukowicz, JPO97
      !!              Bouillon et al., Ocean Modelling 2009
      !!-------------------------------------------------------------------
      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
      !!
      INTEGER ::   ji, jj   ! dummy loop indices
      INTEGER ::   jter     ! local integers
      CHARACTER (len=50) ::   charout

      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
      REAL(wp) ::   zTauO, zTauE, zCor                                       ! temporary scalars

      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
      REAL(wp) ::   zintb, zintn                                             ! dummy argument
      
      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                          ! temporary array for ice strength
      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors
      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                         ! P/delta at T points
      !
      REAL(wp), POINTER, DIMENSION(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points
      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points
      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses
      
      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! shear
      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12                  ! stress tensor components
      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr           ! check convergence
      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                           ! array used for the calculation of ice surface slope:
                                                                             !   ocean surface (ssh_m) if ice is not embedded
                                                                             !   ice top surface if ice is embedded   
      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays
      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice

      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
      !!-------------------------------------------------------------------

      CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )

#if  defined key_lim2 && ! defined key_lim2_vp
# if defined key_agrif
      USE ice_2, vt_s => hsnm
      USE ice_2, vt_i => hicm
# else
      vt_s => hsnm
      vt_i => hicm
# endif
      at_i(:,:) = 1. - frld(:,:)
#endif
#if defined key_agrif && defined key_lim2 
      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
#endif
      !
      !------------------------------------------------------------------------------!
      ! 0) mask at F points for the ice (on the whole domain, not only k_j1,k_jpj) 
      !------------------------------------------------------------------------------!
      ! ocean/land mask
      DO jj = 1, jpjm1
         DO ji = 1, jpim1      ! NO vector opt.
            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
         END DO
      END DO
      CALL lbc_lnk( zfmask, 'F', 1._wp )

      ! Lateral boundary conditions on velocity (modify zfmask)
      zwf(:,:) = zfmask(:,:)
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1   ! vector opt.
            IF( zfmask(ji,jj) == 0._wp ) THEN
               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
            ENDIF
         END DO
      END DO
      DO jj = 2, jpjm1
         IF( zfmask(1,jj) == 0._wp ) THEN
            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
         ENDIF
         IF( zfmask(jpi,jj) == 0._wp ) THEN
            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
         ENDIF
      END DO
      DO ji = 2, jpim1
         IF( zfmask(ji,1) == 0._wp ) THEN
            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
         ENDIF
         IF( zfmask(ji,jpj) == 0._wp ) THEN
            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
         ENDIF
      END DO
      CALL lbc_lnk( zfmask, 'F', 1._wp )

      !------------------------------------------------------------------------------!
      ! 1) define some variables and initialize arrays
      !------------------------------------------------------------------------------!
      ! ecc2: square of yield ellipse eccenticrity
      ecc2    = rn_ecc * rn_ecc
      z1_ecc2 = 1._wp / ecc2

      ! Time step for subcycling
      zdtevp   = rdt_ice / REAL( nn_nevp )
      z1_dtevp = 1._wp / zdtevp

      ! alpha parameters (Bouillon 2009)
#if defined key_lim3
      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
#else
      zalph1 = ( 2._wp * telast ) * z1_dtevp
#endif
      zalph2 = zalph1 * z1_ecc2

      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
      z1_alph2 = 1._wp / ( zalph2 + 1._wp )

      ! Initialise stress tensor 
      zs1 (:,:) = stress1_i (:,:) 
      zs2 (:,:) = stress2_i (:,:)
      zs12(:,:) = stress12_i(:,:)

      ! Ice strength
#if defined key_lim3
      CALL lim_itd_me_icestrength( nn_icestr )
      zpresh(:,:) = tmask(:,:,1) *  strength(:,:)
#else
      zpresh(:,:) = tmask(:,:,1) *  pstar * vt_i(:,:) * EXP( -c_rhg * (1. - at_i(:,:) ) )
#endif

      ! scale factors
      DO jj = k_j1+1, k_jpj-1
         DO ji = fs_2, fs_jpim1
            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
         END DO
      END DO
            
      !
      !------------------------------------------------------------------------------!
      ! 2) Wind / ocean stress, mass terms, coriolis terms
      !------------------------------------------------------------------------------!

      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
         !                                            
         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
         !
         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
         !
         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
         !
      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
         zpice(:,:) = ssh_m(:,:)
      ENDIF

      DO jj = k_j1+1, k_jpj-1
         DO ji = fs_2, fs_jpim1

            ! ice fraction at U-V points
            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)

            ! Ice/snow mass at U-V points
            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
            zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
            zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)

            ! Ocean currents at U-V points
            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
            
            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)

            ! Coriolis at T points (m*f)
            zmf(ji,jj)      = zm1 * fcor(ji,jj)

            ! m/dt
            zmU_t(ji,jj)    = zmassU * z1_dtevp
            zmV_t(ji,jj)    = zmassV * z1_dtevp

            ! Drag ice-atm.
            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)

            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)

            ! masks
            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice

            ! switches
            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin

         END DO
      END DO
      CALL lbc_lnk( zmf, 'T', 1. )
      !
      !------------------------------------------------------------------------------!
      ! 3) Solution of the momentum equation, iterative procedure
      !------------------------------------------------------------------------------!
      !
      !                                               !----------------------!
      DO jter = 1 , nn_nevp                           !    loop over jter    !
         !                                            !----------------------!        
         IF(ln_ctl) THEN   ! Convergence test
            DO jj = k_j1, k_jpj-1
               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
               zv_ice(:,jj) = v_ice(:,jj)
            END DO
         ENDIF

         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
         DO jj = k_j1, k_jpj-1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
            DO ji = 1, jpim1

               ! shear at F points
               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)   &
                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
                  &         ) * r1_e12f(ji,jj) * zfmask(ji,jj)

            END DO
         END DO
         CALL lbc_lnk( zds, 'F', 1. )

         DO jj = k_j1+1, k_jpj-1
            DO ji = 2, jpim1 ! no vector loop

               ! shear**2 at T points (doc eq. A16)
               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
                  &   ) * 0.25_wp * r1_e12t(ji,jj)
              
               ! divergence at T points
               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
                  &    ) * r1_e12t(ji,jj)
               zdiv2 = zdiv * zdiv
               
               ! tension at T points
               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)   &
                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
                  &   ) * r1_e12t(ji,jj)
               zdt2 = zdt * zdt
               
               ! delta at T points
               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )  

               ! P/delta at T points
               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl )
               
               ! stress at T points
               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
             
            END DO
         END DO
         CALL lbc_lnk( zp_delt, 'T', 1. )

         DO jj = k_j1, k_jpj-1
            DO ji = 1, jpim1

               ! P/delta at F points
               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) )
               
               ! stress at F points
               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2

            END DO
         END DO
         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
 
         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
         DO jj = k_j1+1, k_jpj-1
            DO ji = fs_2, fs_jpim1               

               ! U points
               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
                  &                    ) * r1_e2u(ji,jj)                                                                      &
                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
                  &                  ) * r1_e12u(ji,jj)

               ! V points
               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
                  &                    ) * r1_e1v(ji,jj)                                                                      &
                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
                  &                  ) * r1_e12v(ji,jj)

               ! u_ice at V point
               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
               
               ! v_ice at U point
               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)

            END DO
         END DO
         !
         ! --- Computation of ice velocity --- !
         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
         !  Bouillon et al. 2009 (eq 34-35) => stable
         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
            
            DO jj = k_j1+1, k_jpj-1
               DO ji = fs_2, fs_jpim1

                  ! tau_io/(v_oce - v_ice)
                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )

                  ! Coriolis at V-points (energy conserving formulation)
                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )

                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
                  
                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part)
                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
                     &           ) * zmaskV(ji,jj)
               END DO
            END DO
            CALL lbc_lnk( v_ice, 'V', -1. )
            
#if defined key_agrif && defined key_lim2
            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
#endif
#if defined key_bdy
            CALL bdy_ice_lim_dyn( 'V' )
#endif         

            DO jj = k_j1+1, k_jpj-1
               DO ji = fs_2, fs_jpim1
                               
                  ! tau_io/(u_oce - u_ice)
                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )

                  ! Coriolis at U-points (energy conserving formulation)
                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
                  
                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )

                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part)
                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
                     &           ) * zmaskU(ji,jj)
               END DO
            END DO
            CALL lbc_lnk( u_ice, 'U', -1. )
            
#if defined key_agrif && defined key_lim2
            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
#endif
#if defined key_bdy
            CALL bdy_ice_lim_dyn( 'U' )
#endif         

         ELSE ! odd iterations

            DO jj = k_j1+1, k_jpj-1
               DO ji = fs_2, fs_jpim1
                               
                  ! tau_io/(u_oce - u_ice)
                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )

                  ! Coriolis at U-points (energy conserving formulation)
                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  &
                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
                  
                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )

                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part)
                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
                     &           ) * zmaskU(ji,jj)
               END DO
            END DO
            CALL lbc_lnk( u_ice, 'U', -1. )
            
#if defined key_agrif && defined key_lim2
            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
#endif
#if defined key_bdy
            CALL bdy_ice_lim_dyn( 'U' )
#endif         

           DO jj = k_j1+1, k_jpj-1
               DO ji = fs_2, fs_jpim1

                  ! tau_io/(v_oce - v_ice)
                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )

                  ! Coriolis at V-points (energy conserving formulation)
                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  &
                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )

                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
                  
                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part)
                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
                     &           ) * zmaskV(ji,jj)
               END DO
            END DO
            CALL lbc_lnk( v_ice, 'V', -1. )
            
#if defined key_agrif && defined key_lim2
            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
#endif
#if defined key_bdy
            CALL bdy_ice_lim_dyn( 'V' )
#endif         

         ENDIF
         
         IF(ln_ctl) THEN   ! Convergence test
            DO jj = k_j1+1, k_jpj-1
               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
            END DO
            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
         ENDIF
         !
         !                                                ! ==================== !
      END DO                                              !  end loop over jter  !
      !                                                   ! ==================== !
      !
      !------------------------------------------------------------------------------!
      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
      !------------------------------------------------------------------------------!
      DO jj = k_j1, k_jpj-1 
         DO ji = 1, jpim1

            ! shear at F points
            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)   &
               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
               &         ) * r1_e12f(ji,jj) * zfmask(ji,jj)

         END DO
      END DO           
      CALL lbc_lnk( zds, 'F', 1. )
      
      DO jj = k_j1+1, k_jpj-1 
         DO ji = 2, jpim1 ! no vector loop
            
            ! tension**2 at T points
            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)   &
               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
               &   ) * r1_e12t(ji,jj)
            zdt2 = zdt * zdt
            
            ! shear**2 at T points (doc eq. A16)
            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  &
               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  &
               &   ) * 0.25_wp * r1_e12t(ji,jj)
            
            ! shear at T points
            shear_i(ji,jj) = SQRT( zdt2 + zds2 )

            ! divergence at T points
            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
               &            ) * r1_e12t(ji,jj)
            
            ! delta at T points
            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )  
            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
            delta_i(ji,jj) = zdelta + rn_creepl * rswitch

         END DO
      END DO
      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
      
      ! --- Store the stress tensor for the next time step --- !
      stress1_i (:,:) = zs1 (:,:)
      stress2_i (:,:) = zs2 (:,:)
      stress12_i(:,:) = zs12(:,:)

      !
      !------------------------------------------------------------------------------!
      ! 5) Control prints of residual and charge ellipse
      !------------------------------------------------------------------------------!
      !
      ! print the residual for convergence
      IF(ln_ctl) THEN
         WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
         CALL prt_ctl_info(charout)
         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
      ENDIF

      ! print charge ellipse
      ! This can be desactivated once the user is sure that the stress state
      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
      IF(ln_ctl) THEN
         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
            CALL prt_ctl_info(charout)
            DO jj = k_j1+1, k_jpj-1
               DO ji = 2, jpim1
                  IF (zpresh(ji,jj) > 1.0) THEN
                     zsig1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
                     zsig2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
                     CALL prt_ctl_info(charout)
                  ENDIF
               END DO
            END DO
            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
            CALL prt_ctl_info(charout)
         ENDIF
      ENDIF
      !
      CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )

   END SUBROUTINE lim_rhg

#else
   !!----------------------------------------------------------------------
   !!   Default option          Dummy module           NO LIM sea-ice model
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
   END SUBROUTINE lim_rhg
#endif

   !!==============================================================================
END MODULE limrhg