Index: NEMO/trunk/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 10510)
+++ NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 10511)
@@ -861,5 +861,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
 
          !                 !---------------------!
@@ -887,5 +887,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
-         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d    (1:npti), wfx_pnd    (:,:) )
          !
       END SELECT
Index: NEMO/trunk/src/ICE/icethd.F90
===================================================================
--- NEMO/trunk/src/ICE/icethd.F90	(revision 10510)
+++ NEMO/trunk/src/ICE/icethd.F90	(revision 10511)
@@ -152,5 +152,5 @@
             ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2)
             zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 
-            qsb_ice_bot(ji,jj) = rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
+            qsb_ice_bot(ji,jj) = rswitch * rau0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
 
             qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) )
Index: NEMO/trunk/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/trunk/src/OCE/SBC/sbcblk.F90	(revision 10510)
+++ NEMO/trunk/src/OCE/SBC/sbcblk.F90	(revision 10511)
@@ -46,5 +46,5 @@
    USE lib_fortran    ! to use key_nosignedzero
 #if defined key_si3
-   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su, rn_cnd_s, hfx_err_dif
+   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
    USE icethd_dh      ! for CALL ice_thd_snwblow
 #endif
@@ -1099,5 +1099,5 @@
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
-      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat
+      REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat
       !
       ! ECHAM6 constants
@@ -1127,4 +1127,9 @@
       !!----------------------------------------------------------------------
 
+      ! mean temperature
+      WHERE( at_i_b(:,:) > 1.e-20 )   ;   ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
+      ELSEWHERE                       ;   ztm_su(:,:) = rt0
+      ENDWHERE
+      
       ! Momentum Neutral Transfert Coefficients (should be a constant)
       zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
@@ -1137,14 +1142,14 @@
      
       ! Atmospheric and Surface Variables
-      zst(:,:)     = sst_m(:,:) + rt0                                       ! convert SST from Celcius to Kelvin
-      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)  , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
-      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
+      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
+      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
+      zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
       !
       DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
          DO ji = fs_2, fs_jpim1
             ! Virtual potential temperature [K]
-            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
-            zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
-            zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
+            zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
+            zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
+            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
             
             ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
