Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8498)
@@ -332,6 +332,4 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)          [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_thd     !: ice concentration tendency (thermodynamics) [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics)       [s-1]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
@@ -479,22 +477,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmi_dyn  !: Change in ice mass due to ice dynamics (kg/m2/s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dms_dyn  !: Change in snow mass due to ice dynamics (kg/m2/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xmtrp_ice !: X-component of ice mass transport (kg/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_ymtrp_ice !: Y-component of ice mass transport (kg/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xmtrp_snw !: X-component of snow mass transport (kg/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_ymtrp_snw !: Y-component of snow mass transport (kg/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xatrp    !: X-component of area transport (m2/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_yatrp    !: Y-component of area transport (m2/s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_bo    !: Bottom conduction flux (W/m2)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_su    !: Surface conduction flux (W/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_utau_oi  !: X-direction ocean-ice stress
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vtau_oi  !: Y-direction ocean-ice stress  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dssh_dx  !: X-direction sea-surface tilt term (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dssh_dy  !: X-direction sea-surface tilt term (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_corstrx  !: X-direction coriolis stress (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_corstry  !: Y-direction coriolis stress (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_intstrx  !: X-direction internal stress (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_intstry  !: Y-direction internal stress (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sig1     !: Average normal stress in sea ice   
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sig2     !: Maximum shear stress in sea ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_shear    !: Maximum shear of sea-ice velocity field
 
@@ -534,5 +516,5 @@
          &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
          &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
-         &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) , rn_amax_2d(jpi,jpj),   &
+         &      afx_tot(jpi,jpj) , rn_amax_2d(jpi,jpj),                                         &
          &      fhtur  (jpi,jpj) , qlead  (jpi,jpj) ,                                           &
          &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  &
@@ -616,13 +598,5 @@
       ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    & 
                 diag_dmi_dyn(jpi,jpj) , diag_dms_dyn(jpi,jpj) ,    &
-                diag_xmtrp_ice(jpi,jpj), diag_ymtrp_ice(jpi,jpj),  &
-                diag_xmtrp_snw(jpi,jpj), diag_ymtrp_snw(jpi,jpj),  &
-                diag_xatrp(jpi,jpj)    , diag_yatrp(jpi,jpj)    ,  &
                 diag_fc_bo(jpi,jpj)   , diag_fc_su(jpi,jpj)   ,    &
-                diag_utau_oi(jpi,jpj) , diag_vtau_oi(jpi,jpj) ,    &
-                diag_dssh_dx(jpi,jpj) , diag_dssh_dy(jpi,jpj) ,    &
-                diag_corstrx(jpi,jpj) , diag_corstry(jpi,jpj) ,    &
-                diag_intstrx(jpi,jpj) , diag_intstry(jpi,jpj) ,    &
-                diag_sig1(jpi,jpj)    , diag_sig2(jpi,jpj)    ,    &
                 diag_shear(jpi,jpj)   , STAT = ierr(ii) )
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8498)
@@ -27,11 +27,10 @@
    REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
 
-!!gm real variable name starting with a "c" NOT DOCTOR !!!!
    INTEGER  ::   albd_init = 0       ! control flag for initialization  
-   REAL(wp) , PARAMETER ::   c1        = 0.05    ! snow thickness (only for nn_ice_alb=0)
-   REAL(wp) , PARAMETER ::   c2        = 0.10    !  "        "
-   REAL(wp) , PARAMETER ::   rcloud    = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
-   REAL(wp) , PARAMETER ::   r1_c1 = 1. / c1
-   REAL(wp) , PARAMETER ::   r1_c2 = 1. / c2
+   REAL(wp) , PARAMETER ::   rc1    = 0.05    ! snow thickness (only for nn_ice_alb=0)
+   REAL(wp) , PARAMETER ::   rc2    = 0.10    !  "        "
+   REAL(wp) , PARAMETER ::   rcloud = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
+   REAL(wp) , PARAMETER ::   r1_c1 = 1. / rc1
+   REAL(wp) , PARAMETER ::   r1_c2 = 1. / rc2
  
    !                             !!* namelist namsbc_alb *
@@ -114,8 +113,7 @@
       IF ( .NOT. ld_pnd ) THEN   !--- No melt ponds OR radiatively inactive melt ponds
          ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth
-         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = rn_alb_imlt
-                                                       ! !!! MV I think we could replace rt0_ice by rt0 and get rid of rt0
-         ELSE WHERE                                              ;   zalb(:,:,:) = rn_alb_idry
-         END  WHERE
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   ;   zalb(:,:,:) = rn_alb_imlt
+         ELSEWHERE                                           ;   zalb(:,:,:) = rn_alb_idry
+         END WHERE
       ENDIF
 
@@ -127,10 +125,10 @@
          !
          !                       ! Thickness-dependent bare ice albedo
-         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
-         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
-         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
-            &                                                                 - 0.8608 * ph_ice * ph_ice     &
-            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
-         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
+         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_it = zalb
+         ELSEWHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
+         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )   ;   zalb_it = 0.2467 + 0.7049 * ph_ice              &
+            &                                                                  - 0.8608 * ph_ice * ph_ice     &
+            &                                                                  + 0.3812 * ph_ice * ph_ice * ph_ice
+         ELSEWHERE                                        ;   zalb_it = 0.1    + 3.6    * ph_ice
          END WHERE
          !
@@ -139,5 +137,5 @@
             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd ) 
             !                          ! Snow-free ice albedo is a function of pond fraction
-            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   
+            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   
                zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd 
             END WHERE
@@ -150,6 +148,6 @@
                DO ji = 1, jpi
                   !                    ! Freezing snow
-                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
-                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
+                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > rc2
+                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - rc1 ) ) )
                   zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
                      &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * r1_c1  )   &
@@ -157,6 +155,6 @@
                      !
                   !                    ! Melting snow
-                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
+                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > rc2
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - rc2 ) )
                   zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 )   &
                      &      +         zswitch   *   rn_alb_smlt 
@@ -211,5 +209,5 @@
             !
             !                    ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions
-            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )
+            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )
                zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd
             END WHERE
@@ -256,11 +254,11 @@
          !                       !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time)
          zafrac_pnd = 0._wp
-         IF ( ld_pnd ) THEN  
-            WHERE( ph_snw == 0._wp ) ;  zafrac_pnd = pafrac_pnd ;  END WHERE  ! Snow fully "shades" melt ponds
+         IF ( ld_pnd ) THEN
+            WHERE( ph_snw == 0._wp )   zafrac_pnd = pafrac_pnd   ! Snow fully "shades" melt ponds
          ENDIF         
          !
          !                       !--- Specific snow fraction (for now, prescribed)
-         WHERE     ( ph_snw > 0._wp     ) ;  zafrac_snw = 1.
-         ELSE WHERE                       ;  zafrac_snw = 0.
+         WHERE( ph_snw > 0._wp )   ;   zafrac_snw = 1.
+         ELSEWHERE                 ;   zafrac_snw = 0.
          END WHERE
          !
@@ -273,8 +271,8 @@
          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
          z1_c2 = 1. / 0.05
-         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_ice = zalb
-         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_ice = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
+         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_ice = zalb
+         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_ice = zalb       + ( 0.18 - zalb       ) * z1_c1 *  &
             &                                                                       ( LOG(1.5) - LOG(ph_ice) )
-         ELSE WHERE                                       ;  zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
+         ELSEWHERE                                        ;   zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice
          END WHERE
          !
@@ -282,6 +280,6 @@
          z1_c2 = 1. / 0.03
          !
-         WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 );
-         ELSE WHERE                 ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 );
+         WHERE( pt_ice < rt0_snow )   ;   zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 )
+         ELSEWHERE                    ;   zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 )
          END WHERE
          !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8498)
@@ -26,4 +26,5 @@
    USE lib_mpp        ! MPP library
    USE timing         ! Timing
+   USE iom            !
 
    IMPLICIT NONE
@@ -53,4 +54,5 @@
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       REAL(wp) ::   zsal, zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc
+      REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag
       !!----------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('icecor')
@@ -65,69 +67,32 @@
       !
       !
-      !                             !-----------------------------------------------------
-      IF( kn == 2 ) THEN            !  thickness of the smallest category above himin    !
+      IF( kn == 2 ) THEN
          !                          !-----------------------------------------------------
-         !
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-!!gm replace this
-               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes
-               ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch
-!!gm by more readable coding (not slower coding since already a IF in the loop):
-!               IF( a_i(ji,jj,1) >= epsi20 )   ht_i(ji,jj,1) = v_i (ji,jj,1) / a_i(ji,jj,1)
-!!gm
-               IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN
-                  a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin
-               ENDIF
-            END DO
-         END DO
+         !                          !  thickness of the smallest category above himin    !
+         !                          !-----------------------------------------------------
+         WHERE( a_i(:,:,1) >= epsi20 )   ;   ht_i(:,:,1) = v_i (:,:,1) / a_i(:,:,1)
+         ELSEWHERE                       ;   ht_i(:,:,1) = 0._wp
+         END WHERE
+         WHERE( ht_i(:,:,1) < rn_himin )     a_i (:,:,1) = a_i (:,:,1) * ht_i(:,:,1) / rn_himin
          !
       ENDIF
-
       !                             !-----------------------------------------------------
-      at_i(:,:) = a_i(:,:,1)        !  ice concentration should not exceed amax          !
-      DO jl = 2, jpl                !-----------------------------------------------------
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
+      !                             !  ice concentration should not exceed amax          !
+      !                             !-----------------------------------------------------
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
+      DO jl  = 1, jpl
+         WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:)
       END DO
-      !
-!!gm   Question   it seams to me that we have the following equality (dropping the "(ji,jj)":
-!      ( 1. - ( 1. - rn_amax_2d / at_i ) ) =  ( 1. - ( at_i - rn_amax_2d ) / at_i )
-!                                          =  ( at_i - ( at_i - rn_amax_2d ) ) / at_i
-!                                          =  ( + rn_amax_2d  ) / at_i
-!                                          =  rn_amax_2d / at_i
-!     No ?  if yes see "!!gm   better" juste below 
-!gm
-      DO jl  = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
-                  a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
-!!gm   better:    a_i(ji,jj,jl) = a_i(ji,jj,jl) * rn_amax_2d(ji,jj) / at_i(ji,jj)           
-               ENDIF
-            END DO
-         END DO
-      END DO
-!!gm  Other question:  why testing a_i(ji,jj,jl) > 0._wp ?   a_i is >=0, a multiplication by 0 does not change the results....
-!!gm  so at the end, the loop can be recoded without IF as:
-!      WHERE( at_i(:,:) > rn_amax_2d(:,:) )
-!         DO jl  = 1, jpl
-!            a_i(:,:,jl) = a_i(:,:,jl) * MAX( rn_amax_2d(:,:), at_i(:,:) ) / at_i(:,:)
-!         END DO
-!      END WHERE
-!!gm  No?
     
       !                             !-----------------------------------------------------
-      IF (  nn_icesal == 2  ) THEN  !  Ice salinity bounds                               !
+      IF (  nn_icesal == 2  ) THEN  !  salinity stays in bounds [Simin,Simax]            !
       !                             !-----------------------------------------------------
          zzc = rhoic * r1_rdtice
-         DO jl = 1, jpl                  ! salinity stays in bounds [Simin,Simax]
+         DO jl = 1, jpl
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  IF( v_i(ji,jj,jl) > 0._wp ) THEN   ! clem: useless IF ???
-                     zsal = smv_i(ji,jj,jl)
-                     smv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , smv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  )
-                     ! associated salt flux
-                     sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc
-                  ENDIF
+                  zsal = smv_i(ji,jj,jl)
+                  smv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , smv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  )
+                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux
                END DO
             END DO
@@ -157,8 +122,4 @@
          END DO
          CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions
-         !
-!!gm  I think masking here is unnecessary, u_ice already masked and we only introduce zeros in the field
-         u_ice(:,:) = u_ice(:,:) * umask(:,:,1)                            ! mask velocities
-         v_ice(:,:) = v_ice(:,:) * vmask(:,:,1)
       ENDIF
 
@@ -172,7 +133,4 @@
       !                             !-----------------------------------------------------
       CASE( 1 )                        !--- dyn trend diagnostics
-         DO jl  = 1, jpl
-            afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice
-         END DO
          !
 !!gm   here I think the number of ice cat is too small to use a SUM instruction...
@@ -188,12 +146,12 @@
             END DO
          END DO
+         !                       ! concentration tendency (dynamics)
+         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
+         afx_tot(:,:) = zafx(:,:)
+         IF( iom_use('afxdyn') )   CALL iom_put( 'afxdyn' , zafx(:,:) )
          !
       CASE( 2 )                        !--- thermo trend diagnostics & ice aging
          !
-         DO jl  = 1, jpl
-            oa_i(:,:,jl) = oa_i(:,:,jl) +   a_i(:,:,jl)                   * rdt_ice       ! ice natural aging incrementation
-            afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice     ! thermo tendancy
-         END DO
-         afx_tot(:,:) = afx_thd(:,:) + afx_dyn(:,:)
+         oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice   ! ice natural aging incrementation
          !
 !!gm   here I think the number of ice cat is too small to use a SUM instruction...
@@ -209,4 +167,9 @@
             END DO
          END DO
+         !                       ! concentration tendency (total + thermo)
+         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice
+         afx_tot(:,:) = afx_tot(:,:) + zafx(:,:)
+         IF( iom_use('afxthd') )   CALL iom_put( 'afxthd' , zafx(:,:) )
+         IF( iom_use('afxtot') )   CALL iom_put( 'afxtot' , afx_tot(:,:) )
          !
       END SELECT
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8498)
@@ -71,9 +71,5 @@
       REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
       !!----------------------------------------------------------------------
-
-!!gm  Note that glo_sum for a 2D or 3D array use a multiplication by tmask_i(ji,jj)
-!!    so below  the  * tmask(:,:,1) is useless   ===>> I have removed them
-!!    I also move the conversion factor from then glo_sum argument (become a single multiplication
- 
+      !
       IF( icount == 0 ) THEN
          !                          ! salt flux
@@ -271,6 +267,4 @@
                !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj) 
                !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj)
-               !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj)
-               !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj)
                !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj)
                !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj)
@@ -556,6 +550,4 @@
                WRITE(numout,*) ' utau       : ', utau    (ji,jj) 
                WRITE(numout,*) ' vtau       : ', vtau    (ji,jj)
-               WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj)
-               WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj)
             ENDIF
             
@@ -669,5 +661,4 @@
       CALL prt_ctl(tab2d_1=utau       , clinfo1= ' utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
       CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
-      CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
       
    END SUBROUTINE ice_prt3D
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8498)
@@ -65,5 +65,5 @@
       IF( kt == nit000 .AND. lwp ) THEN
          WRITE(numout,*)
-         WRITE(numout,*)'ice_forcing_tau'
+         WRITE(numout,*)'ice_forcing_tau : Surface boundary condition for sea ice (momentum)'
          WRITE(numout,*)'~~~~~~~~~~~~~~~'
       ENDIF
@@ -124,5 +124,5 @@
       IF( kt == nit000 .AND. lwp ) THEN
          WRITE(numout,*)
-         WRITE(numout,*)'ice_forcing_flx'
+         WRITE(numout,*)'ice_forcing_flx : Surface boundary condition for sea ice (flux)'
          WRITE(numout,*)'~~~~~~~~~~~~~~~'
       ENDIF
@@ -132,13 +132,6 @@
 
       ! albedo depends on cloud fraction because of non-linear spectral effects
-      DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
 !!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument !
-               alb_ice(ji,jj,jl) = ( 1. - cldf_ice ) * zalb_cs(ji,jj,jl) + cldf_ice * zalb_os(ji,jj,jl)
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( alb_ice, 'T', 1. )
+      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
       
       SELECT CASE( ksbc )      !==  fluxes over sea ice  ==!
@@ -198,52 +191,52 @@
       INTEGER  ::   jl      ! dummy loop index
       !
-      REAL(wp), DIMENSION(jpi,jpj) :: zalb_m    ! Mean albedo over all categories
-      REAL(wp), DIMENSION(jpi,jpj) :: ztem_m    ! Mean temperature over all categories
-      !
-      REAL(wp), DIMENSION(jpi,jpj) :: z_qsr_m   ! Mean solar heat flux over all categories
-      REAL(wp), DIMENSION(jpi,jpj) :: z_qns_m   ! Mean non solar heat flux over all categories
-      REAL(wp), DIMENSION(jpi,jpj) :: z_evap_m  ! Mean sublimation over all categories
-      REAL(wp), DIMENSION(jpi,jpj) :: z_dqn_m   ! Mean d(qns)/dT over all categories
-      REAL(wp), DIMENSION(jpi,jpj) :: z_devap_m ! Mean d(evap)/dT over all categories
+      REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   ! inverse of concentration
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qsr_m   ! Mean solar heat flux over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qns_m   ! Mean non solar heat flux over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_evap_m  ! Mean sublimation over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_dqn_m   ! Mean d(qns)/dT over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_devap_m ! Mean d(evap)/dT over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zalb_m    ! Mean albedo over all categories
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztem_m    ! Mean temperature over all categories
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('ice_flx_dist')
       !
+      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:)
+      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp
+      END WHERE
+      
       SELECT CASE( k_limflx )       !==  averaged on all ice categories  ==!
       !
       CASE( 0 , 1 )
-         z_qns_m  (:,:) = fice_ice_ave( pqns_ice  (:,:,:) )
-         z_qsr_m  (:,:) = fice_ice_ave( pqsr_ice  (:,:,:) )
-         z_dqn_m  (:,:) = fice_ice_ave( pdqn_ice  (:,:,:) )
-         z_evap_m (:,:) = fice_ice_ave( pevap_ice (:,:,:) )
-         z_devap_m(:,:) = fice_ice_ave( pdevap_ice(:,:,:) )
-!!gm faster coding
-!    REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   ! 
-! ...
-!      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:)
-!      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp
-!      END WHERE
-!      z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
-!      z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
-!      z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
-!      z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
-!      z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
-!! and remove the 2 functions: fice_ice_ave and fice_cell_ave
-!!gm
+         !
+         ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) )  
+         !
+         z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
+         z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
          DO jl = 1, jpl
-            pdqn_ice  (:,:,jl) = z_dqn_m  (:,:)
-            pdevap_ice(:,:,jl) = z_devap_m(:,:)
             pqns_ice  (:,:,jl) = z_qns_m (:,:)
             pqsr_ice  (:,:,jl) = z_qsr_m (:,:)
+            pdqn_ice  (:,:,jl) = z_dqn_m  (:,:)
             pevap_ice (:,:,jl) = z_evap_m(:,:)
+            pdevap_ice(:,:,jl) = z_devap_m(:,:)
          END DO
          !
+         DEALLOCATE( z_qns_m, z_qsr_m, z_dqn_m, z_evap_m, z_devap_m )  
+         !
       END SELECT
       !
       SELECT CASE( k_limflx )       !==  redistribution on all ice categories  ==!
+      !
       CASE( 1 , 2 )
          !
-         zalb_m(:,:) = fice_ice_ave( palb_ice(:,:,:) )
-         ztem_m(:,:) = fice_ice_ave( ptn_ice (:,:,:) )
+         ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) )  
+         !
+         zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:)
+         ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:)
          DO jl = 1, jpl
             pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
@@ -252,4 +245,6 @@
          END DO
          !
+         DEALLOCATE( zalb_m, ztem_m )  
+         !
       END SELECT
       !
@@ -257,34 +252,4 @@
       !
    END SUBROUTINE ice_flx_dist
-
-!!gm TO BE REMOVED ====>>>>>
-   FUNCTION fice_cell_ave ( ptab )
-      !!--------------------------------------------------------------------------
-      !! * Compute average over categories, for grid cell (ice covered and free ocean)
-      !!--------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) :: ptab
-      REAL(wp), DIMENSION(jpi,jpj)                 :: fice_cell_ave
-      INTEGER :: jl
-      !!--------------------------------------------------------------------------
-      fice_cell_ave(:,:) = a_i(:,:,1) * ptab (:,:,1)
-      DO jl = 2, jpl
-         fice_cell_ave(:,:) = fice_cell_ave(:,:) + a_i(:,:,jl) * ptab (:,:,jl)
-      END DO
-   END FUNCTION fice_cell_ave
-
-
-   FUNCTION fice_ice_ave ( ptab )
-      !!--------------------------------------------------------------------------
-      !! * Compute average over categories, for ice covered part of grid cell
-      !!--------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) ::   ptab   !
-      REAL(wp), DIMENSION(jpi,jpj)                 :: fice_ice_ave
-      !!--------------------------------------------------------------------------
-      WHERE ( at_i (:,:) > 0.0_wp ) ; fice_ice_ave (:,:) = fice_cell_ave( ptab (:,:,:) ) / at_i (:,:)
-      ELSEWHERE                     ; fice_ice_ave (:,:) = 0.0_wp
-      END WHERE
-   END FUNCTION fice_ice_ave
-
-!!gm <<<<<<====  end of TO BE REMOVED 
 
 #else
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90	(revision 8498)
@@ -546,18 +546,19 @@
          WRITE(numout,*)
          WRITE(numout,*) 'ice_ist_init : ice parameters inititialisation '
-         WRITE(numout,*) '~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_limini     = ', ln_limini
-         WRITE(numout,*) '   ice initialization from a netcdf file      ln_limini_file  = ', ln_limini_file
-         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
-         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
-         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
-         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
-         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
-         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
-         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
-         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
-         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
-         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
-         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namiceini'
+         WRITE(numout,*) '      initialization with ice (T) or not (F)       ln_limini     = ', ln_limini
+         WRITE(numout,*) '      ice initialization from a netcdf file      ln_limini_file  = ', ln_limini_file
+         WRITE(numout,*) '      threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
+         WRITE(numout,*) '      initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
+         WRITE(numout,*) '      initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
+         WRITE(numout,*) '      initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
+         WRITE(numout,*) '      initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
+         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
+         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
+         WRITE(numout,*) '      initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
+         WRITE(numout,*) '      initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
+         WRITE(numout,*) '      initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
+         WRITE(numout,*) '      initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
       ENDIF
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8498)
@@ -426,11 +426,4 @@
             ENDIF
             !
-!!gm use of rswitch not faster as there is already IF in the DO-loop ==>>> use IF which is more readable
-!            rswitch    = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) )
-!            zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch
-!            !
-!            rswitch    = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) )
-!            zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch        
-
             IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
             ELSE                                  ;   zworkv(ji) = 0._wp
@@ -439,5 +432,4 @@
             ELSE                                  ;   zworka(ji) = 0._wp
             ENDIF
-!!gm end
             !
             a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8498)
@@ -48,14 +48,9 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
-
+   !
    REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
    REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer
-
-!!gm Cp is  1) not DOCTOR, 
-!!          2) misleading name: heat capacity instead of a constant,
-!!          3) recomputed at each time-step, whereas it is stored in the module memory...
-!!             ===>>> compute it one for all inside the IF( kt == nit000 ) (i.e. without the ".AND. lwp")
-   REAL(wp)            ::   Cp                  ! ??? !!gm  Not doctor ! 
-   
+   !
+   REAL(wp)            ::   zdrho               ! 
    !
    !
@@ -111,6 +106,5 @@
       INTEGER  ::   niter              ! local integer 
       INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging
-      REAL(wp) ::   za, zfac, zcs_2    ! local scalar
-      CHARACTER (len = 15) ::   fieldid
+      REAL(wp) ::   za                 ! local scalar
       REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s)
       !                                                 ! (ridging ice area - area of new ridges) / dt
@@ -125,22 +119,12 @@
       IF( nn_timing == 1 )  CALL timing_start('icerdgrft')
 
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*)'icerdgrft : ice ridging and rafting'
-         WRITE(numout,*)'~~~~~~~~~~'
-      ENDIF
-!!gm should be:      
-!      IF( kt == nit000 ) THEN
-!         IF(lwp) WRITE(numout,*)
-!         IF(lwp) WRITE(numout,*)'icerdgrft : ???'
-!         IF(lwp) WRITE(numout,*)'~~~~~~~~~~'
-!         !
-!         Cp    = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
-!         !
-!!gm why not adding here zcs_2 computation
-!         !
-!      ENDIF
-!!gm end
-      
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)'icerdgrft : ice ridging and rafting'
+         IF(lwp) WRITE(numout,*)'~~~~~~~~~~'
+         !
+         zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
+         !
+      ENDIF      
       !                    ! conservation test
       IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
@@ -149,6 +133,4 @@
       ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
       !-----------------------------------------------------------------------------!
-      Cp    = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE
-      zcs_2 = rn_cs * 0.5_wp
       !
       CALL ice_rdgrft_ridgeprep                             ! prepare ridging
@@ -175,5 +157,5 @@
             !  (thick, newly ridged ice).
 
-            closing_net(ji,jj) = zcs_2 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
+            closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
 
             ! 2.2 divu_adv
@@ -220,8 +202,6 @@
                za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
                IF    ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i
-                  zfac          = - ato_i(ji,jj) / za
                   opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
                ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
-                  zfac          = ( asum(ji,jj) - ato_i(ji,jj) ) / za
                   opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
                ENDIF
@@ -238,6 +218,5 @@
                   za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
                   IF( za  >  a_i(ji,jj,jl) ) THEN
-                     zfac = a_i(ji,jj,jl) / za
-                     closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
+                     closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za
                   ENDIF
                END DO
@@ -254,5 +233,5 @@
          asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
 
-         ! 3.5 Do we keep on iterating ???
+         ! 3.5 Do we keep on iterating?
          !-----------------------------------------------------------------------------!
          ! Check whether asum = 1.  If not (because the closing and opening
@@ -316,34 +295,17 @@
       !!---------------------------------------------------------------------!
       INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar     !!gm DOCTOR norme should start with z !!!!!
-      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
+      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zdummy   ! local scalar
+      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   zGsum     ! zGsum(n) = sum of areas in categories 0 to n
       !------------------------------------------------------------------------------!
 
-      Gstari        = 1._wp / rn_gstar    
-      astari        = 1._wp / rn_astar    
-      aksum(:,:)    = 0._wp
-      athorn(:,:,:) = 0._wp
-      aridge(:,:,:) = 0._wp
-      araft (:,:,:) = 0._wp
+      z1_gstar      = 1._wp / rn_gstar    
+      z1_astar      = 1._wp / rn_astar    
 
       CALL ice_var_zapsmall   ! Zero out categories with very small areas
 
-!      DO jl = 1, jpl          ! Ice thickness needed for rafting
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!!gm               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
-!!gm               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-!               IF( a_i(ji,jj,jl) >= epsi20 ) THEN   ;   ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl)
-!               ELSE                                 ;   ht_i(ji,jj,jl) = 0._wp
-!               ENDIF
-!           END DO
-!         END DO
-!      END DO
-!!gm or even better :
-!     !                       ! Ice thickness needed for rafting
+      !                       ! Ice thickness needed for rafting
       WHERE( a_i(:,:,:) >= epsi20 )   ;   ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:)
       ELSEWHERE                       ;   ht_i(:,:,:) = 0._wp
       END WHERE
-!!gm end
 
       !------------------------------------------------------------------------------!
@@ -356,17 +318,11 @@
       !
       ! Compute cumulative thickness distribution function
-      ! Compute the cumulative thickness distribution function Gsum,
-      ! where Gsum(n) is the fractional area in categories 0 to n.
+      ! Compute the cumulative thickness distribution function zGsum,
+      ! where zGsum(n) is the fractional area in categories 0 to n.
       ! initial value (in h = 0) equals open water area
-      Gsum(:,:,-1) = 0._wp
-      Gsum(:,:,0 ) = ato_i(:,:)
-      ! for each value of h, you have to add ice concentration then
+      zGsum(:,:,-1) = 0._wp
+      zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:)
       DO jl = 1, jpl
-         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
-      END DO
-      !
-      ! Normalize the cumulative distribution to 1
-      DO jl = 0, jpl
-         Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:)
+         zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl)
       END DO
 
@@ -388,10 +344,10 @@
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  IF    ( Gsum(ji,jj,jl)   < rn_gstar ) THEN
-                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
-                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
-                  ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN
-                     athorn(ji,jj,jl) = Gstari * ( rn_gstar       - Gsum(ji,jj,jl-1) ) *  &
-                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar       ) * Gstari )
+                  IF    ( zGsum(ji,jj,jl)   < rn_gstar ) THEN
+                     athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * &
+                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar )
+                  ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN
+                     athorn(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  &
+                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar        ) * z1_gstar )
                   ELSE
                      athorn(ji,jj,jl) = 0._wp
@@ -403,10 +359,10 @@
       ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
          !                        
-         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
+         zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) )        ! precompute exponential terms using zGsum as a work array
          DO jl = -1, jpl
-            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
+            zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy
          END DO
          DO jl = 0, jpl
-             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
+            athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl)
          END DO
          !
@@ -418,18 +374,18 @@
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  zdummy           = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) )
-                  aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl)
+                  aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl)
                   araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl)
                END DO
             END DO
          END DO
-      ELSEIF( ln_ridging .AND. .NOT.ln_rafting ) THEN   !- ridging alone
-         DO jl = 1, jpl
-            aridge(:,:,jl) = athorn(:,:,jl)
-         END DO
-      ELSEIF( ln_rafting .AND. .NOT.ln_ridging ) THEN   !- rafting alone   
-         DO jl = 1, jpl
-            araft(:,:,jl) = athorn(:,:,jl)
-         END DO
+      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
+         aridge(:,:,:) = athorn(:,:,1:jpl)
+         araft (:,:,:) = 0._wp
+      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
+        aridge(:,:,:) = 0._wp
+         araft (:,:,:) = athorn(:,:,1:jpl)
+      ELSE                                               !- no ridging & no rafting
+         aridge(:,:,:) = 0._wp
+         araft (:,:,:) = 0._wp         
       ENDIF
 
@@ -441,13 +397,13 @@
       ! 
       ! This parameterization is a modified version of Hibler (1980).
-      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
+      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
       !  and for very thick ridging ice must be >= krdgmin*hi
       !
       ! The minimum ridging thickness, hrmin, is equal to 2*hi 
       !  (i.e., rafting) and for very thick ridging ice is
-      !  constrained by hrmin <= (hrmean + hi)/2.
+      !  constrained by hrmin <= (zhmean + hi)/2.
       ! 
       ! The maximum ridging thickness, hrmax, is determined by
-      !  hrmean and hrmin.
+      !  zhmean and hrmin.
       !
       ! These modifications have the effect of reducing the ice strength
@@ -466,9 +422,9 @@
             DO ji = 1, jpi
                IF ( athorn(ji,jj,jl) > 0._wp ) THEN
-                  hrmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
-                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) )
-                  hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl)
+                  zhmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
+                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) )
+                  hrmax(ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl)
                   hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
-                  krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 )
+                  krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 )
                   !
                   ! Normalization factor : aksum, ensures mass conservation
@@ -498,16 +454,13 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
-      !
-      CHARACTER (len=80) ::   fieldid   ! field identifier     !!gm DOCTOR name wrong
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area retreats, excluding area of new ridges
       !
       INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
       INTEGER ::   ij                ! horizontal index, combines i and j loops
       INTEGER ::   icells            ! number of cells with a_i > puny
-      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration
-      REAL(wp) ::   zwfx_snw         ! snow mass flux increment
+      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration and new area going to jl2
 
       INTEGER , DIMENSION(jpij) ::   indxi, indxj   ! compressed indices
-      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to n2
+      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to jl2
 
       REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged 
@@ -656,8 +609,8 @@
             !  During the next time step, thermo_rates will determine whether
             !  the ocean cools or new ice grows.
-            zwfx_snw           = ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
-               &                 + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
-
-            wfx_snw_dyn(ji,jj)  =   wfx_snw_dyn(ji,jj) + zwfx_snw
+            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
+               &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
+
+!!!            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + zwfx_snw
 
             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
@@ -722,5 +675,5 @@
 
                ! update jl1
-               e_i  (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
+               e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
 
             END DO
@@ -806,6 +759,6 @@
       !
       INTEGER             ::   ji,jj, jl   ! dummy loop indices
-      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation    !!gm not DOCTOR : start with i !!!
-      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing    !!gm not DOCTOR : start with i !!!
+      INTEGER             ::   ismooth     ! smoothing the resistance to deformation
+      INTEGER             ::   itframe     ! number of time steps for the P smoothing
       REAL(wp)            ::   zp, z1_3    ! local scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
@@ -813,68 +766,39 @@
       !!----------------------------------------------------------------------
 
-      !------------------------------------------------------------------------------!
-      ! 1) Initialize
-      !------------------------------------------------------------------------------!
-      strength(:,:) = 0._wp
-
-      !------------------------------------------------------------------------------!
-      ! 2) Compute thickness distribution of ridging and ridged ice
-      !------------------------------------------------------------------------------!
-      CALL ice_rdgrft_ridgeprep
-
-      !------------------------------------------------------------------------------!
-      ! 3) Rothrock(1975)'s method
-      !------------------------------------------------------------------------------!
-      IF( kstrngth == 1 ) THEN
+      !                              !--------------------------------------------------!
+      CALL ice_rdgrft_ridgeprep      ! Thickness distribution of ridging and ridged ice !
+      !                              !--------------------------------------------------!
+
+      !                              !--------------------------------------------------!
+      IF( kstrngth == 1 ) THEN       ! Ice strength => Rothrock (1975) method           !
+      !                              !--------------------------------------------------!
          z1_3 = 1._wp / 3._wp
          DO jl = 1, jpl
-            DO jj= 1, jpj
-               DO ji = 1, jpi
-                  !
-                  IF( athorn(ji,jj,jl) > 0._wp ) THEN
-                     !----------------------------
-                     ! PE loss from deforming ice
-                     !----------------------------
-                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
-
-                     !--------------------------
-                     ! PE gain from rafting ice
-                     !--------------------------
-                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
-
-                     !----------------------------
-                     ! PE gain from ridging ice
-                     !----------------------------
-                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 *  &
-                        &                              ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) +         &
-                        &                                hrmin(ji,jj,jl) * hrmin(ji,jj,jl) +         &
-                        &                                hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )  
-                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                      
-                  ENDIF
-                  !
-               END DO
-            END DO
-         END DO
-   
-         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1)
-                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
-         ksmooth = 1
-
-      !------------------------------------------------------------------------------!
-      ! 4) Hibler (1979)' method
-      !------------------------------------------------------------------------------!
-      ELSE                      ! kstrngth ne 1:  Hibler (1979) form
-         !
+            WHERE( athorn(:,:,jl) > 0._wp )
+               strength(:,:) =  -         athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice
+                  &             + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice
+                  &             +         aridge(:,:,jl) * krdg(:,:,jl) * z1_3 *   &        ! PE gain from ridging ice
+                  &                      ( hrmax(:,:,jl) * hrmax(:,:,jl) +         &
+                  &                        hrmin(:,:,jl) * hrmin(:,:,jl) +         &
+                  &                        hrmax(:,:,jl) * hrmin(:,:,jl) )
+            ELSEWHERE
+               strength(:,:) = 0._wp
+            END WHERE
+         END DO
+         strength(:,:) = rn_pe_rdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)
+                         ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
+         ismooth = 1
+      !                              !--------------------------------------------------!
+      ELSE                           ! Ice strength => Hibler (1979) method             !
+      !                              !--------------------------------------------------!
          strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
          !
-         ksmooth = 1
-         !
-      ENDIF                     ! kstrngth
-      !
-      !------------------------------------------------------------------------------!
-      ! 5) Impact of brine volume
-      !------------------------------------------------------------------------------!
+         ismooth = 1
+         !
+      ENDIF
+      !                              !--------------------------------------------------!
+      IF( ln_icestr_bvf ) THEN       ! Impact of brine volume                           !
+      !                              !--------------------------------------------------!
       ! CAN BE REMOVED
-      IF( ln_icestr_bvf ) THEN
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -883,12 +807,9 @@
          END DO
       ENDIF
-      !
-      !------------------------------------------------------------------------------!
-      ! 6) Smoothing ice strength
-      !------------------------------------------------------------------------------!
-      SELECT CASE( ksmooth )
-      !                       !-------------------
-      CASE( 1 )               ! Spatial smoothing
-         !                    !-------------------
+      !                              !--------------------------------------------------!
+      SELECT CASE( ismooth )         ! Smoothing ice strength                           !
+      !                              !--------------------------------------------------!
+      !
+      CASE( 1 )               !--- Spatial smoothing
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -903,5 +824,5 @@
             END DO
          END DO
-
+         
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -911,7 +832,5 @@
          CALL lbc_lnk( strength, 'T', 1. )
          !
-         !                    !--------------------
-      CASE( 2 )               ! Temporal smoothing
-         !                    !--------------------
+      CASE( 2 )               !--- Temporal smoothing
          IF ( kt_ice == nit000 ) THEN
             zstrp1(:,:) = 0._wp
@@ -922,8 +841,8 @@
             DO ji = 2, jpim1
                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 
-                  numts_rm = 1 ! number of time steps for the running mean
-                  IF ( zstrp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
-                  IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
-                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm
+                  itframe = 1 ! number of time steps for the running mean
+                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
+                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
+                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
                   zstrp2  (ji,jj) = zstrp1  (ji,jj)
                   zstrp1  (ji,jj) = strength(ji,jj)
@@ -932,5 +851,5 @@
             END DO
          END DO
-         CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions
+         CALL lbc_lnk( strength, 'T', 1. )
          !
       END SELECT
@@ -970,20 +889,21 @@
       IF (lwp) THEN                          ! control print
          WRITE(numout,*)
-         WRITE(numout,*)'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
-         WRITE(numout,*)'~~~~~~~~~~~~~~~'
-         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
-         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
-         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
-         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
-         WRITE(numout,*)'   Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
-         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
-         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
-         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
-         WRITE(numout,*)'   Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
-         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
-         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
-         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft  
-         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
-         WRITE(numout,*)'   Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
+         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for mechanical ice redistribution '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namiceitdme'
+         WRITE(numout,*) '      Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
+         WRITE(numout,*) '      Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
+         WRITE(numout,*) '      Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
+         WRITE(numout,*) '      Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
+         WRITE(numout,*) '      Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging
+         WRITE(numout,*) '      Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
+         WRITE(numout,*) '      Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
+         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
+         WRITE(numout,*) '      Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
+         WRITE(numout,*) '      Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
+         WRITE(numout,*) '      Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
+         WRITE(numout,*) '      Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft  
+         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
+         WRITE(numout,*) '      Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
       ENDIF
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8498)
@@ -31,4 +31,5 @@
    USE prtctl         ! Print control
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE iom            !
 #if defined key_agrif
    USE agrif_lim3_interp
@@ -111,5 +112,4 @@
       INTEGER ::   ji, jj       ! dummy loop indices
       INTEGER ::   jter         ! local integers
-      CHARACTER (len=50) ::   charout
 
       REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
@@ -121,8 +121,8 @@
       REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! 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) ::   zfac_x, zfac_y
+      REAL(wp) ::   zshear, zdum1, zdum2
       
       REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
@@ -152,4 +152,23 @@
       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
+      !! --- diags
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zswi, zsig1, zsig2, zsig3
+      !! --- SIMIP diags
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress  
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)      
       !!-------------------------------------------------------------------
 
@@ -671,88 +690,126 @@
 
       !------------------------------------------------------------------------------!
-      ! 5) SIMIP diagnostics
-      !------------------------------------------------------------------------------!
-
-!!gm  encapsulate with a flag (iom_use of the variable or better a flag defined one for all testing if one of the
-!!    diag is output.  NB the diag_...  are should only be ALLOCATED if the flag is true !
-
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-             rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
-
-             ! Stress tensor invariants (normal and shear stress N/m)
-             diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
-             diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
-
-             ! Stress terms of the momentum equation (N/m2)
-             diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
-             diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
-
-             diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
-             diag_corstry(ji,jj) = zCory(ji,jj) * rswitch
-
-             diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
-             diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
-           
-             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
-             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
-
-             ! 2D ice mass, snow mass, area transport arrays (X, Y)
-             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
-             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
-
-             diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
-             diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
-
-             diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
-             diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
-
-             diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
-             diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
-
-         END DO
-      END DO
-
-      CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   &
-         &                  diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   &
-         &                  diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   & 
-         &                  diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    )
-
-      CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    )
-
-      CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1.,   &
-         &                  diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1.,   &
-         &                  diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.    )
-
-      !
-      !------------------------------------------------------------------------------!
-      ! 6) Control prints of residual and charge ellipse
-      !------------------------------------------------------------------------------!
-      !
-      ! print the residual for convergence
-      IF(ln_ctl) THEN
-         WRITE(charout,FMT="('ice_rhg_evp  : res =',D23.16, ' iter =',I4)") zresm, jter
-         CALL prt_ctl_info(charout)
-         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
+      ! 5) diagnostics
+      !------------------------------------------------------------------------------!
+      ! --- ellipse --- !
+      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
          !
-         ! 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. (2008) for more details
-         IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN
-            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. '
-            CALL prt_ctl_info(charout)
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF (strength(ji,jj) > 1.0) THEN
-                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
-                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )
-                     WRITE(charout,FMT="('ice_rhg_evp  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
-                     CALL prt_ctl_info(charout)
-                  ENDIF
-               END DO
-            END DO
-            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. '
-            CALL prt_ctl_info(charout)
-         ENDIF
+         ALLOCATE( zswi(jpi,jpj) , zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
+         !
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+            END DO
+         END DO
+         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
+                  &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  &
+                  &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
+
+               zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
+
+               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
+
+!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
+!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
+!!               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
+!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
+               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
+               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                    ! shear stress
+               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
+         !
+         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
+         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
+         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
+         !
+         DEALLOCATE( zswi , zsig1 , zsig2 , zsig3 )
+      ENDIF
+      
+      ! --- SIMIP --- !
+      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
+         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
+         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
+         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
+
+         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
+            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
+            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
+            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
+         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+               
+               ! Stress tensor invariants (normal and shear stress N/m)
+               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
+               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
+               
+               ! Stress terms of the momentum equation (N/m2)
+               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
+               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
+               
+               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
+               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
+               
+               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
+               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
+               
+               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
+               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
+               
+               ! 2D ice mass, snow mass, area transport arrays (X, Y)
+               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
+               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
+               
+               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
+               zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
+               
+               zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
+               zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
+               
+               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
+               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
+               
+            END DO
+         END DO
+         
+         CALL lbc_lnk_multi(   zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
+            &                  zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
+            &                  zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
+            &                  zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
+         
+         CALL lbc_lnk_multi(   zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1.    )
+         
+         CALL lbc_lnk_multi(   zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
+            &                  zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
+            &                  zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
+         
+         IF ( iom_use( 'normstr'  ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
+         IF ( iom_use( 'sheastr'  ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
+         IF ( iom_use( 'dssh_dx'  ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
+         IF ( iom_use( 'dssh_dy'  ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
+         IF ( iom_use( 'corstrx'  ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
+         IF ( iom_use( 'corstry'  ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
+         IF ( iom_use( 'intstrx'  ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
+         IF ( iom_use( 'intstry'  ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
+         IF ( iom_use( 'utau_oi'  ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
+         IF ( iom_use( 'vtau_oi'  ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
+         IF ( iom_use( 'xmtrpice' ) )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
+         IF ( iom_use( 'ymtrpice' ) )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport 
+         IF ( iom_use( 'xmtrpsnw' ) )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
+         IF ( iom_use( 'ymtrpsnw' ) )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
+         IF ( iom_use( 'xatrp'    ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
+         IF ( iom_use( 'yatrp'    ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
+
+         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
+            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
+            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
+            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
+
       ENDIF
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8498)
@@ -122,20 +122,16 @@
          !                                   !-----------------------!
          !
-         kt_ice = kt       ! Ice model time step
+         kt_ice = kt                              ! -- Ice model time step
          !
 # if defined key_agrif
          IF( .NOT. Agrif_Root() )  lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
 # endif
-
-         !                 ! mean surface ocean current at ice velocity point
-         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
-         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
-!!gm the umask, vmask above is useless as ssu_m, ssv_m are build from masked un,vn (used t be here for B-grid sea-ice)
-
-         !                 ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
-         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
+         u_oce(:,:) = ssu_m(:,:)                  ! -- mean surface ocean current
+         v_oce(:,:) = ssv_m(:,:)
+         !
+         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )   ! -- freezing temperature [Kelvin] (set to rt0 over land)
          t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
          !
-         CALL ice_bef      ! Store previous ice values
+         CALL store_fields                        ! -- Store now ice values
          !
          !------------------------------------------------!
@@ -179,5 +175,5 @@
          CALL ice_var_glo2eqv            ! ht_i and ht_s for ice albedo calculation
          CALL ice_var_agg(1)             ! at_i for coupling 
-         CALL ice_bef                    ! Store previous ice values
+         CALL store_fields               ! Store now ice values
 
          !------------------------------------------------------!
@@ -223,14 +219,14 @@
 !!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
 !!# endif
-         IF( ln_limdiahsb )           CALL ice_dia( kt )        ! -- Diagnostics and outputs 
-         !
-                                      CALL ice_wri( 1 )         ! -- Ice outputs 
+         IF( ln_limdiahsb )         CALL ice_dia( kt )        ! -- Diagnostics and outputs 
+         !
+                                    CALL ice_wri( 1 )         ! -- Ice outputs 
          !
          IF( kt == nit000 .AND. ln_rstart )   &                !!gm  I don't understand the ln_rstart, if needed, add a comment, please
-            &                         CALL iom_close( numrir )  ! close input ice restart file
-         !
-         IF( lrst_ice )               CALL ice_rst_write( kt )  ! -- Ice restart file 
-         !
-         IF( ln_limctl )              CALL ice_ctl( kt )        ! alerts in case of model crash
+            &                       CALL iom_close( numrir )  ! close input ice restart file
+         !
+         IF( lrst_ice )             CALL ice_rst_write( kt )  ! -- Ice restart file 
+         !
+         IF( ln_limctl )            CALL ice_ctl( kt )        ! alerts in case of model crash
          !
       ENDIF   ! End sea-ice time step only
@@ -241,5 +237,5 @@
       ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
       !    using before instantaneous surf. currents
-      IF( ln_limdyn )                 CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )
+      IF( ln_limdyn )               CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
       !
@@ -476,7 +472,7 @@
 
 
-   SUBROUTINE ice_bef
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_bef  ***
+   SUBROUTINE store_fields
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE store_fields  ***
       !!
       !! ** purpose :  store ice variables at "before" time step
@@ -525,5 +521,5 @@
       CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. )
       !
-   END SUBROUTINE ice_bef
+   END SUBROUTINE store_fields
 
 
@@ -570,5 +566,4 @@
             !
             afx_tot(ji,jj) = 0._wp   ;
-            afx_dyn(ji,jj) = 0._wp   ;   afx_thd(ji,jj) = 0._wp
             !
             diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8498)
@@ -193,14 +193,9 @@
       !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema) 
       !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            hfx_out(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean               
-               &             - qlead(ji,jj) * r1_rdtice                                    &  ! heat flux taken from the ocean where there is open water ice formation
-               &             - at_i(ji,jj) * fhtur(ji,jj)                                  &  ! heat flux taken by turbulence
-               &             - at_i(ji,jj) *  fhld(ji,jj)                                     ! heat flux taken during bottom growth/melt 
-                                                                                              !    (fhld should be 0 while bott growth)
-         END DO
-      END DO
-
+      hfx_out(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean               
+         &           - qlead(:,:) * r1_rdtice                                &  ! heat flux taken from the ocean where there is open water ice formation
+         &           - at_i (:,:) * fhtur(:,:)                               &  ! heat flux taken by turbulence
+         &           - at_i (:,:) *  fhld(:,:)                                  ! heat flux taken during bottom growth/melt 
+                                                                                !    (fhld should be 0 while bott growth)
       !-------------------------------------------------------------------------------------------!
       ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
@@ -288,20 +283,19 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ji, jk   ! dummy loop indices
-      REAL(wp) ::   ztmelts, zaaa, zbbb, zccc, zdiscrim  ! local scalar 
+      REAL(wp) ::   ztmelts, z1_2cp, zbbb, zccc  ! local scalar 
       !!-------------------------------------------------------------------
       ! Recover ice temperature
+      z1_2cp = 1._wp / ( 2._wp * cpic )
       DO jk = 1, nlay_i
          DO ji = 1, nidx
-            ztmelts       =  -tmut * s_i_1d(ji,jk) + rt0
+            ztmelts       = -tmut * s_i_1d(ji,jk)
             ! Conversion q(S,T) -> T (second order equation)
-            zaaa          =  cpic
-            zbbb          =  ( rcp - cpic ) * ( ztmelts - rt0 ) + e_i_1d(ji,jk) * r1_rhoic - lfus
-            zccc          =  lfus * ( ztmelts - rt0 )
-            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
-            t_i_1d(ji,jk) =  rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
+            zbbb          = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus
+            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) )
+            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * z1_2cp
             
             ! mask temperature
-            rswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
-            t_i_1d(ji,jk) =  rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
+            rswitch       = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
+            t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0
          END DO 
       END DO 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90	(revision 8498)
@@ -32,6 +32,4 @@
 CONTAINS
 
-!!gm  even comment line of more than 130 character may cause compilation errors
-!!gm         ===>>> reformat the text
    SUBROUTINE ice_thd_da
       !!-------------------------------------------------------------------
@@ -46,6 +44,8 @@
       !!                   W = m1 * (Tw -Tf)**m2                    --- originally from Josberger 1979 ---
       !!                      (Tw - Tf) = elevation of water temp above freezing
-      !!                      m1 and m2 = (1.6e-6 , 1.36) best fit from field experiment near the coast of Prince Patrick Island (Perovich 1983) => static ice
-      !!                      m1 and m2 = (3.0e-6 , 1.36) best fit from MIZEX 84 experiment (Maykut and Perovich 1987) => moving ice
+      !!                      m1 and m2 = (1.6e-6 , 1.36) best fit from field experiment near the coast of Prince Patrick Island
+      !!                                                                                           (Perovich 1983) => static ice
+      !!                      m1 and m2 = (3.0e-6 , 1.36) best fit from MIZEX 84 experiment
+      !!                                                                                (Maykut and Perovich 1987) => moving ice
       !!
       !!                   P = N * pi * D                           --- from Rothrock and Thorndike 1984 ---
@@ -59,5 +59,6 @@
       !!                      Astar = 1 / ( 1 - (Dmin/Dmax)**(1/beta) )
       !!                      Dmin = minimum floe diameter (recommended to be 8m +- 20%)
-      !!                      Dmax = maximum floe diameter (recommended to be 300m, but it does not impact melting much except for Dmax<100m)
+      !!                      Dmax = maximum floe diameter (recommended to be 300m,
+      !!                                                    but it does not impact melting much except for Dmax<100m)
       !!                      beta = 1.0 +-20% (recommended value)
       !!                           = 0.3 best fit for western Fram Strait and Antarctica
@@ -113,7 +114,9 @@
          ! --- Calculate reduction of total sea ice concentration --- !
          zdfloe = rn_dmin * ( zastar / ( zastar - at_i_1d(ji) ) )**rn_beta         ! Mean floe caliper diameter [m]
-         zperi  = at_i_1d(ji) * rpi / ( zcs * zdfloe )                             ! Mean perimeter of the floe = N*pi*D = (A/cs*D^2)*pi*D [m.m-2]
+         !
+         zperi  = at_i_1d(ji) * rpi / ( zcs * zdfloe )                             ! Mean perimeter of the floe [m.m-2]
+         !                                                                         !    = N*pi*D = (A/cs*D^2)*pi*D
          zwlat  = zm1 * ( MAX( 0._wp, sst_1d(ji) - ( t_bo_1d(ji) - rt0 ) ) )**zm2  ! Melt speed rate [m/s]
-         
+         !
          zda_tot(ji) = MIN( zwlat * zperi * rdt_ice, at_i_1d(ji) )                 ! sea ice concentration decrease (>0)
       
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8498)
@@ -102,5 +102,5 @@
       REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing
 
-      REAL(wp) :: zswitch_sal
+      REAL(wp) ::   zswitch_sal
 
       INTEGER  ::   num_iter_max      ! Heat conservation 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8498)
@@ -53,14 +53,9 @@
       !!---------------------------------------------------------------------
 
-      !--------------------------------------------------------------------|
-      ! 1) salinity constant in time                                       |
-      !--------------------------------------------------------------------|
-      ! do nothing
-
-      !----------------------------------------------------------------------|
-      !  2) salinity varying in time                                         |
-      !----------------------------------------------------------------------|
-      IF(  nn_icesal == 2  ) THEN
-
+      SELECT CASE ( nn_icesal )
+      !
+      !              !---------------------------------------------!
+      CASE( 2 )      !  time varying salinity with linear profile  !
+      !              !---------------------------------------------!
          DO ji = 1, nidx
 
@@ -96,11 +91,11 @@
          CALL ice_var_salprof1d
          !
-      ENDIF 
-
-      !------------------------------------------------------------------------------|
-      !  3) vertical profile of salinity, constant in time                           |
-      !------------------------------------------------------------------------------|
-      IF(  nn_icesal == 3  )   CALL ice_var_salprof1d
+      !                 !---------------------------------------------!
+      CASE( 3 )         ! constant salinity with a fix profile        ! (Schwarzacher (1959) multiyear salinity profile(mean = 2.30)
+      !                 !---------------------------------------------!
+         CALL ice_var_salprof1d
       !
+   END SELECT
+   !
    END SUBROUTINE ice_thd_sal
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8498)
@@ -86,8 +86,8 @@
       !!  
       !! ** Purpose :   Update the surface ocean boundary condition for heat 
-      !!              salt and mass over areas where sea-ice is non-zero
+      !!                salt and mass over areas where sea-ice is non-zero
       !!         
       !! ** Action  : - computes the heat and freshwater/salt fluxes
-      !!              at the ice-ocean interface.
+      !!                at the ice-ocean interface.
       !!              - Update the ocean sbc
       !!     
@@ -133,15 +133,7 @@
          DO ji = 1, jpi
 
-            !------------------------------------------!
-            !      heat flux at the ocean surface      !
-            !------------------------------------------!
             ! Solar heat flux reaching the ocean = zqsr (W.m-2) 
             !---------------------------------------------------
-            zqsr = qsr_tot(ji,jj)
-            DO jl = 1, jpl
-               zqsr = zqsr - a_i_b(ji,jj,jl) * (  qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) 
-            END DO
-!!gm  why not like almost everywhere else :
-!!gm        zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * (  qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) )
+            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) )
 
             ! Total heat flux reaching the ocean = hfx_out (W.m-2) 
@@ -204,13 +196,11 @@
       END DO
 
-      !-----------------------------------------------!
-      !   Storing the transmitted variables           !
-      !-----------------------------------------------!
+      ! Storing the transmitted variables
+      !----------------------------------
       fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction            
       tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                      
 
-      !------------------------------------------------------------------------!
-      !    Snow/ice albedo (only if sent to coupler, useless in forced mode)   !
-      !------------------------------------------------------------------------!
+      ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
+      !------------------------------------------------------------------
       CALL ice_alb( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8498)
@@ -157,8 +157,8 @@
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   ze_i, z1_CpR, zdiscrim, zaaa, z1_2aaa             ! local scalars
-      REAL(wp) ::   ze_s, zL_Cp , ztmelts , zbbb, zccc                !   -      -
-      REAL(wp) ::   zhmax, z1_zhmax, zsm_i, zcpMcpic, zt_i   !   -      -
-      REAL(wp) ::   zlay_i, zlay_s   !   -      -
+      REAL(wp) ::   ze_i, z1_cp, z1_2cp             ! local scalars
+      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
+      REAL(wp) ::   zhmax, z1_zhmax, zsm_i          !   -      -
+      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
       !!------------------------------------------------------------------
@@ -181,5 +181,5 @@
       ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness
 
-	  zhmax    =          hi_max(jpl)
+      zhmax    =          hi_max(jpl)
       z1_zhmax =  1._wp / hi_max(jpl)               
       WHERE( ht_i(:,:,jpl) > zhmax )               !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area
@@ -205,7 +205,5 @@
       !-------------------
       zlay_i   = REAL( nlay_i , wp )    ! number of layers
-      zaaa     = cpic                   ! Conversion q(S,T) -> T (second order equation)
-      z1_2aaa  = 1._wp / ( 2._wp * zaaa )
-      zcpMcpic = rcp - cpic
+      z1_2cp  = 1._wp / ( 2._wp * cpic )
       DO jl = 1, jpl
          DO jk = 1, nlay_i
@@ -214,15 +212,13 @@
                   IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
                      !
-                     ze_i    =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
-                     ztmelts = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
-                     !
-                     zbbb     =  zcpMcpic * ztmelts + ze_i * r1_rhoic - lfus
-                     zccc     =  lfus * ztmelts
-                     zdiscrim =  SQRT(  MAX( zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp)  )
-                     zt_i     = - ( zbbb + zdiscrim ) * z1_2aaa                          ! [C]
-                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( zt_i , ztmelts )  ) + rt0   ! [K] with bounds: -100 < zt_i < ztmelts
+                     ze_i             =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
+                     ztmelts          = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
+                     ! Conversion q(S,T) -> T (second order equation)
+                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
+                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
+                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * z1_2cp , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
                      !
                   ELSE                                !--- no ice
-                     t_i(ji,jj,jk,jl) =  rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C)
+                     t_i(ji,jj,jk,jl) = rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C)
 !!clem: I think we should impose rt0 instead
                   ENDIF
@@ -235,25 +231,13 @@
       ! Snow temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere)
       !--------------------
-      z1_CpR = 1._wp / ( cpic * rhosn )
-      zL_Cp  = lfus  /   cpic
       zlay_s = REAL( nlay_s , wp )
+      z1_cp  = 1._wp / cpic
       DO jk = 1, nlay_s
          WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
-            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0       
+            t_s(:,:,jk,:) = MAX( -100._wp , MIN( z1_cp * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
          ELSEWHERE                           !--- no ice
             t_s(:,:,jk,:) =  rt0 - 100._wp         ! impose 173.15 K (i.e. -100 C)
          END WHERE
       END DO
-!!gm perhaps better like this (?) :  (if this compile, yes! one test instead of nlay_s tests)
-!      WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
-!         DO jk = 1, nlay_s
-!            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0       
-!         END DO
-!      ELSEWHERE                           !--- no ice
-!         DO jk = 1, nlay_s
-!            t_s(:,:,jk,:) =  rt0 - 100._wp      ! impose 173.15 K (i.e. -100 C)
-!         END DO
-!      END WHERE
-!!gm end better ?
 
       ! integrated values 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8497)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8498)
@@ -54,9 +54,7 @@
       INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2, zmiss_val
-      REAL(wp) ::  zs12, zshear
       REAL(wp), DIMENSION(jpi,jpj)     ::  z2d, zswi, zmiss !  2D workspace
       REAL(wp), DIMENSION(jpi,jpj)     ::  zfb              ! ice freeboard
       REAL(wp), DIMENSION(jpi,jpj)     ::  zamask, zamask15 ! 15% concentration mask
-      REAL(wp), DIMENSION(jpi,jpj)     ::  zsig1, zsig2, zsig3
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::  zswi2, zmiss2
       !
@@ -196,8 +194,4 @@
       CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean      
  
-      CALL iom_put( "afxtot"     , afx_tot              )        ! concentration tendency (total)
-      CALL iom_put( "afxdyn"     , afx_dyn              )        ! concentration tendency (dynamics)
-      CALL iom_put( "afxthd"     , afx_thd              )        ! concentration tendency (thermo)
-
       CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   !  
       CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   !  
@@ -219,34 +213,4 @@
       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
 
-!!gm ====>>>>>  THIS should be moved in icerhg_evp    (generalize this everywhere it is possible and logic...)
-      ! specific outputs for EVP rheology
-      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
-         zsig1(:,:) = 0._wp; zsig2(:,:) = 0._wp; zsig3(:,:) = 0._wp;
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zs12 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
-                  &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  &
-                  &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
-
-               zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zs12 * zs12 ) ! shear stress  
-
-               z2da = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
-
-!!               zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
-!!               zsig2(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
-!!               zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
-!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
-               zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
-               zsig2(ji,jj) = 0.5_wp * z2da * ( zshear )                    ! shear stress
-               zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
-            END DO
-         END DO
-         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
-         CALL iom_put( "isig1" , zsig1 )
-         CALL iom_put( "isig2" , zsig2 )
-         CALL iom_put( "isig3" , zsig3 )
-      ENDIF
-!!gm  <<<<<<======= end
-
 !!gm ====>>>>>  THIS should be moved where at_ip, vt_ip are computed fro the last time in the time-step  (limmpd)
       ! MV MP 2016
@@ -339,32 +303,9 @@
       IF ( iom_use( "vice_mv"  ) ) CALL iom_put( "vice_mv"     ,   v_ice(:,:)               * zswi(:,:) + zmiss(:,:) )   ! ice velocity v component
       
-      IF ( iom_use( "xmtrpice" ) ) CALL iom_put( "xmtrpice"     ,  diag_xmtrp_ice(:,:)      )                            ! X-component of sea-ice mass transport (kg/s)
-      IF ( iom_use( "ymtrpice" ) ) CALL iom_put( "ymtrpice"     ,  diag_ymtrp_ice(:,:)      )                            ! Y-component of sea-ice mass transport 
-
-      IF ( iom_use( "xmtrpsnw" ) ) CALL iom_put( "xmtrpsnw"     ,  diag_xmtrp_snw(:,:)      )                            ! X-component of snow mass transport (kg/s)
-      IF ( iom_use( "ymtrpsnw" ) ) CALL iom_put( "ymtrpsnw"     ,  diag_ymtrp_snw(:,:)      )                            ! Y-component of snow mass transport
-
-      IF ( iom_use( "xatrp"    ) ) CALL iom_put( "xatrp"        ,  diag_xatrp(:,:)              )                        ! X-component of ice area transport
-      IF ( iom_use( "yatrp"    ) ) CALL iom_put( "yatrp"        ,  diag_yatrp(:,:)              )                        ! Y-component of ice area transport
-
       IF ( iom_use( "utau_ice" ) ) CALL iom_put( "utau_ice"     ,  utau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (x)
       IF ( iom_use( "vtau_ice" ) ) CALL iom_put( "vtau_ice"     ,  vtau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (y)
 
-      IF ( iom_use( "utau_oi"  ) ) CALL iom_put( "utau_oi"     ,   diag_utau_oi(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Ocean stress term in force balance (x)
-      IF ( iom_use( "vtau_oi"  ) ) CALL iom_put( "vtau_oi"     ,   diag_vtau_oi(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Ocean stress term in force balance (y)
 
       IF ( iom_use( "icestr"   ) ) CALL iom_put( "icestr"      ,   strength(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Ice strength
-
-      IF ( iom_use( "dssh_dx"  ) ) CALL iom_put( "dssh_dx"     ,   diag_dssh_dx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Sea-surface tilt term in force balance (x)
-      IF ( iom_use( "dssh_dy"  ) ) CALL iom_put( "dssh_dy"     ,   diag_dssh_dy(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Sea-surface tilt term in force balance (y)
-
-      IF ( iom_use( "corstrx"  ) ) CALL iom_put( "corstrx"     ,   diag_corstrx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Coriolis force term in force balance (x)
-      IF ( iom_use( "corstry"  ) ) CALL iom_put( "corstry"     ,   diag_corstry(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Coriolis force term in force balance (y)
-
-      IF ( iom_use( "intstrx"  ) ) CALL iom_put( "intstrx"     ,   diag_intstrx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Internal force term in force balance (x)
-      IF ( iom_use( "intstry"  ) ) CALL iom_put( "intstry"     ,   diag_intstry(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Internal force term in force balance (y)
-
-      IF ( iom_use( "normstr"  ) ) CALL iom_put( "normstr"     ,   diag_sig1(:,:)           * zswi(:,:) + zmiss(:,:) )   ! Normal stress
-      IF ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   diag_sig2(:,:)           * zswi(:,:) + zmiss(:,:) )   ! Shear stress
 
       !--------------------------------
