Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8325)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8326)
@@ -23,5 +23,5 @@
    USE dom_oce        ! ocean space and time domain variables
    USE ice            ! sea-ice variables
-   USE sbc_oce , ONLY : sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl
+   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl
    USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
       &                 fr1_i0, fr2_i0, nn_limflx
@@ -55,5 +55,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -458,4 +458,5 @@
          CALL tab_2d_1d( nbpb, wfx_snw_sub_1d(1:nbpb), wfx_snw_sub  , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, wfx_ice_sub_1d(1:nbpb), wfx_ice_sub  , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, wfx_err_sub_1d(1:nbpb), wfx_err_sub  , jpi, jpj, npb(1:nbpb) )
          !
          CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) )
@@ -487,8 +488,12 @@
          CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, hfx_out_1d (1:nbpb), hfx_out         , jpi, jpj, npb(1:nbpb) )
          !
          ! SIMIP diagnostics
          CALL tab_2d_1d( nbpb, diag_fc_bo_1d   (1:nbpb), diag_fc_bo  , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, diag_fc_su_1d   (1:nbpb), diag_fc_su  , jpi, jpj, npb(1:nbpb) )
+         ! ocean surface fields
+         CALL tab_2d_1d( nbpb, sst_1d(1:nbpb), sst_m, jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, sss_1d(1:nbpb), sss_m, jpi, jpj, npb(1:nbpb) )
          !
       CASE( 2 )            ! from 1D to 2D
@@ -516,4 +521,5 @@
          CALL tab_1d_2d( nbpb, wfx_snw_sub   , npb, wfx_snw_sub_1d(1:nbpb), jpi, jpj )
          CALL tab_1d_2d( nbpb, wfx_ice_sub   , npb, wfx_ice_sub_1d(1:nbpb), jpi, jpj )
+         CALL tab_1d_2d( nbpb, wfx_err_sub   , npb, wfx_err_sub_1d(1:nbpb), jpi, jpj )
          !
          CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj )
@@ -545,4 +551,5 @@
          CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
          CALL tab_1d_2d( nbpb, hfx_err_dif   , npb, hfx_err_dif_1d(1:nbpb), jpi, jpj )
+         CALL tab_1d_2d( nbpb, hfx_out       , npb, hfx_out_1d(1:nbpb)   , jpi, jpj )
          !
          CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8325)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8326)
@@ -18,5 +18,4 @@
    USE par_oce        ! ocean parameters
    USE phycst         ! physical constants (OCE directory) 
-   USE sbc_oce , ONLY : sst_m, sss_m
    USE ice            ! LIM variables
    USE thd_ice        ! LIM thermodynamics
@@ -70,5 +69,4 @@
       !! 
       INTEGER  ::   ji , jk        ! dummy loop indices
-      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
       INTEGER  ::   iter
 
@@ -87,5 +85,4 @@
       REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
       REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
-      REAL(wp) ::   zsstK        ! SST (K)
 
       REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3)
@@ -403,6 +400,5 @@
       ! remaining "potential" evap is sent to ocean
       DO ji = kideb, kiut
-         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
-         wfx_err_sub(ii,ij) = wfx_err_sub(ii,ij) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
+         wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
       END DO
 
@@ -445,7 +441,5 @@
                   &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )
 
-               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
-
-               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
+               s_i_new(ji)        = zswitch_sal * zfracs * sss_1d(ji)  &  ! New ice salinity
                                   + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
                ! New ice growth
@@ -611,7 +605,6 @@
          wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
          !    
-         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
-         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
+         hfx_out_1d(ji)  = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
 
          IF( ln_limctl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
@@ -637,8 +630,6 @@
 
          ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
-         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
          zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0
-         zsstK          = sst_m(ii,ij) + rt0                                
-         zEw            = rcp * ( zsstK - rt0 )
+         zEw            = rcp * sst_1d(ji)
          zQm            = zfmdt * zEw 
          
@@ -647,9 +638,9 @@
 
          ! Contribution to salt flux
-         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 
+         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 
 
          ! virtual salt flux to keep salinity constant
          IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
-            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_m(ii,ij) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean
+            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean
                &                            - sm_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean 
          ENDIF
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 8325)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 8326)
@@ -16,5 +16,4 @@
    USE par_oce        ! ocean parameters
    USE phycst         ! physical constants (ocean directory)
-   USE sbc_oce , ONLY : sss_m
    USE ice            ! LIM variables
    USE thd_ice        ! LIM thermodynamics
@@ -51,5 +50,5 @@
       INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
       !
-      INTEGER  ::   ii, ij, ji, jk               ! dummy loop indices 
+      INTEGER  ::   ji, jk                       ! dummy loop indices 
       REAL(wp) ::   iflush, igravdr              ! local scalars
       REAL(wp) ::   zs_sni, zsm_i_gd, zsm_i_fl, zsm_i_si, zsm_i_bg   ! local scalars
@@ -68,9 +67,8 @@
          DO ji = kideb, kiut
 
-            ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
             !---------------------------------------------------------
             !  Update ice salinity from snow-ice and bottom growth
             !---------------------------------------------------------
-            zs_sni   = sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic   ! Salinity of snow ice
+            zs_sni   = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic   ! Salinity of snow ice
             rswitch  = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
             zsm_i_si = ( zs_sni      - sm_i_1d(ji) ) *             dh_snowice(ji)  / MAX( ht_i_1d(ji), epsi20 ) * rswitch ! snow-ice    
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8325)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8326)
@@ -48,4 +48,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_dif_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_out_1d
 
    ! heat flux associated with ice-atmosphere mass exchange
@@ -62,4 +63,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sub_1d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_ice_sub_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_err_sub_1d 
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_bog_1d    
@@ -121,4 +123,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   diag_fc_su_1d      !: <==> the 2D  diag_fc_su
 
+   ! surface fields from the ocean
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sst_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sss_1d
+
    INTEGER , PUBLIC ::   jiindex_1d   ! 1D index of debugging point
 
@@ -135,5 +141,5 @@
       !!---------------------------------------------------------------------!
       INTEGER ::   thd_ice_alloc   ! return value
-      INTEGER ::   ierr(5), ii
+      INTEGER ::   ierr(6), ii
       !!---------------------------------------------------------------------!
       ierr(:) = 0
@@ -149,5 +155,5 @@
          &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      &
          &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) ,   &
-         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , STAT=ierr(ii) )
+         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , hfx_out_1d(jpij), STAT=ierr(ii) )
       !
       ii = ii + 1
@@ -156,5 +162,5 @@
          &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) ,  &
          &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  &
-         &      wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij) ,                                    &
+         &      wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,              &
          &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         &
          &      qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0         (jpij) ,                     &  
@@ -177,4 +183,7 @@
       ii = ii + 1
       ALLOCATE( diag_fc_bo_1d(jpij)      , diag_fc_su_1d(jpij)      , STAT=ierr(ii) )
+      !
+      ii = ii + 1
+      ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) )
 
       thd_ice_alloc = MAXVAL( ierr(:) )
