Index: /NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/field_def_nemo-ice.xml	(revision 10311)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/field_def_nemo-ice.xml	(revision 10312)
@@ -29,5 +29,6 @@
           <field id="icemask15"    long_name="Ice mask (0 if ice conc. lower than 15%, 1 otherwise)"   standard_name="sea_ice_mask15"                            unit="" />
 	  <field id="icepres"      long_name="Fraction of time steps with sea ice"                     standard_name="sea_ice_time_fraction"                     unit="" />
-	  
+          <field id="fasticepres"  long_name="Fraction of time steps with landfast ice"                standard_name="fast_ice_time_fraction"                    unit="" />
+ 
 	  <!-- general fields -->
           <field id="icemass"      long_name="Sea-ice mass per area"                                   standard_name="sea_ice_amount"                            unit="kg/m2"/>
Index: /NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/namelist_ice_ref	(revision 10311)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/cfgs/SHARED/namelist_ice_ref	(revision 10312)
@@ -52,14 +52,16 @@
    ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection)
    ln_dynADV        = .false.         !  dyn.: only advection w prescribed vel.(rn_uvice + advection)
-      rn_uice       =   0.00001       !        prescribed ice u-velocity
-      rn_vice       =   0.            !        prescribed ice v-velocity
-   rn_ishlat        =   2.            !  free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2)
-   ln_landfast      = .false.         !  landfast ice parameterization (T or F)                           
-      rn_gamma      =   0.15          !     fraction of ocean depth that ice must reach to initiate landfast
-                                      !        recommended range: [0.1 ; 0.25]
-      rn_icebfr     =  10.            !     maximum bottom stress per unit area of contact [N/m2]                 
-                                      !        a very large value ensures ice velocity=0 even with a small contact area
-                                      !        recommended range: ?? (should be greater than atm-ice stress => >0.1 N/m2)
-      rn_lfrelax    =   1.e-5         !     relaxation time scale to reach static friction [s-1]
+      rn_uice       =   0.5           !        prescribed ice u-velocity
+      rn_vice       =   0.5           !        prescribed ice v-velocity
+   rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2)
+   ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016
+   ln_landfast_home = .false.         !  landfast: parameterization from "home made"
+      rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast
+                                      !          recommended range: [0.1 ; 0.25] - L16=0.125 - home=0.15
+      rn_icebfr     =  15.            !        ln_landfast_L16:  maximum bottom stress per unit volume [N/m3]
+                                      !        ln_landfast_home: maximum bottom stress per unit area of contact [N/m2]
+                                      !          recommended range: ?? L16=15 - home=10
+      rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1]
+      rn_tensile    =   0.2           !        isotropic tensile strength
 /
 !------------------------------------------------------------------------------
Index: /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/ice.F90	(revision 10311)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/ice.F90	(revision 10312)
@@ -129,10 +129,15 @@
    !                                     !!** ice-dynamics namelist (namdyn) **
    REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
-   LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F) 
-   REAL(wp), PUBLIC ::   rn_gamma         !:    fraction of ocean depth that ice must reach to initiate landfast ice
-   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (landfast ice) 
-   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice) 
-   !
-   !                                     !!** ice-rheology namelist (namrhg) **
+   LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016 
+   LOGICAL , PUBLIC ::   ln_landfast_home !: landfast ice parameterizationfrom home made 
+   REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice
+   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home) 
+   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction
+   REAL(wp), PUBLIC ::   rn_tensile       !:    isotropic tensile strength
+   !
+   !                                     !!** ice-ridging/rafting namelist (namdyn_rdgrft) **
+   REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength (also used for landfast param)
+   !
+   !                                     !!** ice-rheology namelist (namdyn_rhg) **
    LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F) 
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
Index: /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90	(revision 10311)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn.F90	(revision 10312)
@@ -84,14 +84,12 @@
          WRITE(numout,*)'~~~~~~~'
       ENDIF
-
       !                      
-      IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction
+      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization
          tau_icebfr(:,:) = 0._wp
          DO jl = 1, jpl
-            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_gamma )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
+            WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
          END DO
-         IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr )  
-      ENDIF
-
+      ENDIF
+      !
       zhmax(:,:,:) = h_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
       DO jl = 1, jpl
@@ -244,6 +242,6 @@
       !!
       NAMELIST/namdyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  &
-         &             rn_ishlat  ,                                            &
-         &             ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax
+         &             rn_ishlat ,                                             &
+         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile
       !!-------------------------------------------------------------------
       !
@@ -261,13 +259,15 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdyn:'
-         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL   = ', ln_dynFULL
-         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV = ', ln_dynRHGADV
-         WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV    = ', ln_dynADV
-         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')'
-         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat    = ', rn_ishlat
-         WRITE(numout,*) '      Landfast: param (T or F)                                ln_landfast  = ', ln_landfast
-         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_gamma     = ', rn_gamma
-         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr    = ', rn_icebfr
-         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax   = ', rn_lfrelax
+         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL      = ', ln_dynFULL
+         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV
+         WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV       = ', ln_dynADV
+         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')'
+         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat
+         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16
+         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home
+         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra
+         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr
+         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax
+         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile
          WRITE(numout,*)
       ENDIF
@@ -289,6 +289,10 @@
       ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
       ENDIF
-      !                                      !--- NO Landfast ice : set to zero once for all
-      IF( .NOT.ln_landfast )   tau_icebfr(:,:) = 0._wp
+      !                                      !--- Landfast ice
+      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp
+      !
+      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN
+         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' )
+      ENDIF
       !
       CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters
Index: /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rdgrft.F90	(revision 10311)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rdgrft.F90	(revision 10312)
@@ -39,15 +39,15 @@
 
    ! Variables shared among ridging subroutines
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   closing_net     ! net rate at which area is removed    (1/s)
-      !                                                 ! (ridging ice area - area of new ridges) / dt
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   opning         ! rate of opening due to divergence/shear
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   closing_gross  ! rate at which area removed, not counting area of new ridges
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmin    ! minimum ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hrmax    ! maximum ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hraft    ! thickness of rafted ice
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aridge   ! participating ice ridging
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   araft    ! participating ice rafting
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_net     ! net rate at which area is removed    (1/s)
+      !                                                               ! (ridging ice area - area of new ridges) / dt
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   opning          ! rate of opening due to divergence/shear
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_gross   ! rate at which area removed, not counting area of new ridges
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   apartf          ! participation function; fraction of ridging/closing associated w/ category n
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmin           ! minimum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmax           ! maximum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hraft           ! thickness of rafted ice
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hi_hrdg         ! thickness of ridging ice / mean ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aridge          ! participating ice ridging
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   araft           ! participating ice rafting
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_i_2d
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ze_s_2d
@@ -59,5 +59,4 @@
    LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79)
    REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
-   REAL(wp) ::   rn_crhg          ! determines changes in ice strength
    REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging            
    LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
Index: /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rhg_evp.F90	(revision 10311)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icedyn_rhg_evp.F90	(revision 10312)
@@ -118,7 +118,9 @@
       REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
       REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
-      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
       REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars
+      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
+      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
       !
       REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
@@ -136,4 +138,5 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param)
       REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
@@ -255,5 +258,9 @@
          END DO
       END DO
-            
+
+      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
+      IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile
+      ELSE                                               ;   zkt = 0._wp
+      ENDIF
       !
       !------------------------------------------------------------------------------!
@@ -273,5 +280,5 @@
          zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
          !
-      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
+      ELSE                               !== non-embedded sea ice: use ocean surface for slope calculation ==!
          zpice(:,:) = ssh_m(:,:)
       ENDIF
@@ -328,4 +335,43 @@
       CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )
       !
+      !                                  !== Landfast ice parameterization ==!
+      !
+      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ! ice thickness at U-V points
+               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)
+               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)
+               ! ice-bottom stress at U points
+               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
+               zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               ! ice-bottom stress at V points
+               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
+               zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
+               ! ice_bottom stress at T points
+               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
+               tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
+            END DO
+         END DO
+         CALL lbc_lnk( tau_icebfr(:,:), 'T', 1. )
+         !
+      ELSEIF( ln_landfast_home ) THEN          !-- Home made
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zTauU_ib(ji,jj) = tau_icebfr(ji,jj)
+               zTauV_ib(ji,jj) = tau_icebfr(ji,jj)
+            END DO
+         END DO
+         !
+      ELSE                                     !-- no landfast
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zTauU_ib(ji,jj) = 0._wp
+               zTauV_ib(ji,jj) = 0._wp
+            END DO
+         END DO
+      ENDIF
+      IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) )
+
       !------------------------------------------------------------------------------!
       ! 3) Solution of the momentum equation, iterative procedure
@@ -391,7 +437,7 @@
                ENDIF
                
-               ! 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
+               ! stress at T points (zkt/=0 if landfast)
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
              
             END DO
@@ -412,6 +458,6 @@
                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
+               ! stress at F points (zkt/=0 if landfast)
+               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
 
             END DO
@@ -461,6 +507,6 @@
                   !
                   !                 !--- tau_bottom/v_ice
-                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
+                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = - zTauV_ib(ji,jj) / zvel
                   !
                   !                 !--- Coriolis at V-points (energy conserving formulation)
@@ -473,5 +519,5 @@
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
@@ -509,6 +555,6 @@
                   !
                   !                 !--- tau_bottom/u_ice
-                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
+                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = - zTauU_ib(ji,jj) / zvel
                   !
                   !                 !--- Coriolis at U-points (energy conserving formulation)
@@ -521,5 +567,5 @@
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
@@ -559,6 +605,6 @@
                   !
                   !                 !--- tau_bottom/u_ice
-                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
+                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = - zTauU_ib(ji,jj) / zvel
                   !
                   !                 !--- Coriolis at U-points (energy conserving formulation)
@@ -571,5 +617,5 @@
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
@@ -607,6 +653,6 @@
                   !
                   !                 !--- tau_bottom/v_ice
-                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
-                  ztauB = - tau_icebfr(ji,jj) / zvel
+                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = - zTauV_ib(ji,jj) / zvel
                   !
                   !                 !--- Coriolis at v-points (energy conserving formulation)
@@ -619,5 +665,5 @@
                   !
                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
                   !
                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
Index: /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icewri.F90	(revision 10311)
+++ /NEMO/branches/2018/dev_r9947_SI3_advection/src/ICE/icewri.F90	(revision 10312)
@@ -50,5 +50,5 @@
       INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::   z2da, z2db, zrho1, zrho2
-      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d !  2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d, zfast !  2D workspace
       REAL(wp), DIMENSION(jpi,jpj)     ::   zmsk00, zmsk05, zmsk15, zmsksn ! O%, 5% and 15% concentration mask and snow mask
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zmsk00l, zmsksnl               ! cat masks
@@ -132,5 +132,6 @@
       IF( iom_use('vtau_ai'  ) )   CALL iom_put( "vtau_ai", vtau_ice * zmsk00     )   ! Wind stress term in force balance (y)
 
-      IF( iom_use('icevel') ) THEN 
+      IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN 
+        ! module of ice velocity
          DO jj = 2 , jpjm1
             DO ji = 2 , jpim1
@@ -141,5 +142,11 @@
          END DO
          CALL lbc_lnk( z2d, 'T', 1. )
-         IF( iom_use('icevel') )   CALL iom_put( "icevel" , z2d                   )   ! ice velocity module
+         IF( iom_use('icevel') )   CALL iom_put( "icevel" , z2d )
+
+        ! record presence of fast ice
+         WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp
+         ELSEWHERE                                                ; zfast(:,:) = 0._wp
+         END WHERE
+         IF( iom_use('fasticepres') )   CALL iom_put( "fasticepres" , zfast )
       ENDIF
 
