Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 4177)
@@ -682,17 +682,7 @@
       ! used to prevent the applied increments taking the temperature below the local freezing point 
 
-#if defined key_cice 
-        fzptnz(:,:,:) = -1.8_wp
-#else 
-        DO jk = 1, jpk
-           DO jj = 1, jpj
-              DO ji = 1, jpk
-                 fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) )                   & 
-                                                  - 2.154996e-4_wp *       tsn(ji,jj,jk,jp_sal)   ) * tsn(ji,jj,jk,jp_sal)  & 
-                                                  - 7.53e-4_wp * fsdepw(ji,jj,jk)       ! (pressure in dbar) 
-              END DO
-           END DO
-        END DO
-#endif 
+      DO jk=1, jpkm1
+         fzptnz (:,:,jk) = tfreez( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) )
+      ENDDO
 
       IF ( ln_asmiau ) THEN
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4177)
@@ -21,4 +21,6 @@
    USE bdy_par         ! (for lk_bdy)
    USE timing          ! preformance summary
+   USE lib_fortran
+   USE sbcrnf
 
    IMPLICIT NONE
@@ -33,4 +35,5 @@
    REAL(dp)                                ::   surf_tot   , vol_tot             !
    REAL(dp)                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
+   REAL(dp)                                ::   frc_wn_t      , frc_wn_s ! global forcing trends
    REAL(dp)                                ::   fact1                            ! conversion factors
    REAL(dp)                                ::   fact21    , fact22               !     -         -
@@ -38,4 +41,5 @@
    REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf      , ssh_ini              !
    REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
+   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini
 
    !! * Substitutions
@@ -67,9 +71,13 @@
       INTEGER    ::   jk                          ! dummy loop indice
       REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
+      REAL(dp)   ::   zdiff_hc1   , zdiff_sc1     ! heat and salt content variations of ssh
       REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
+      REAL(dp)   ::   zerr_hc1    , zerr_sc1      ! Non conservation due to free surface
       REAL(dp)   ::   z1_rau0                     ! local scalars
       REAL(dp)   ::   zdeltat                     !    -     -
       REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
       REAL(dp)   ::   z_frc_trd_v                 !    -     -
+      REAL(dp)   ::   z_wn_trd_t , z_wn_trd_s   !    -     -
+      REAL(dp)   ::   z_ssh_hc , z_ssh_sc   !    -     -
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('dia_hsb')
@@ -79,18 +87,27 @@
       ! ------------------------- !
       z1_rau0 = 1.e0 / rau0
-      z_frc_trd_v = z1_rau0 * SUM( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) )     ! volume fluxes
-      z_frc_trd_t =           SUM( sbc_tsc(:,:,jp_tem) * surf(:,:) )     ! heat fluxes
-      z_frc_trd_s =           SUM( sbc_tsc(:,:,jp_sal) * surf(:,:) )     ! salt fluxes
+      z_frc_trd_v = z1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) )     ! volume fluxes
+      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )     ! heat fluxes
+      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )     ! salt fluxes
+      ! Add runoff heat & salt input
+      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) )
+      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) )
       ! Add penetrative solar radiation
-      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * SUM( qsr     (:,:) * surf(:,:) )
+      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) )
       ! Add geothermal heat flux
-      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * SUM( qgh_trd0(:,:) * surf(:,:) )
-      IF( lk_mpp ) THEN
-         CALL mpp_sum( z_frc_trd_v )
-         CALL mpp_sum( z_frc_trd_t )
-      ENDIF
+      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +  glob_sum( qgh_trd0(:,:) * surf(:,:) )
+      IF( .NOT. lk_vvl ) THEN
+         z_wn_trd_t = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) )
+         z_wn_trd_s = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) )
+      ENDIF
+
       frc_v = frc_v + z_frc_trd_v * rdt
       frc_t = frc_t + z_frc_trd_t * rdt
       frc_s = frc_s + z_frc_trd_s * rdt
+      !                                          ! Advection flux through fixed surface (z=0)
+      IF( .NOT. lk_vvl ) THEN
+         frc_wn_t = frc_wn_t + z_wn_trd_t * rdt
+         frc_wn_s = frc_wn_s + z_wn_trd_s * rdt
+      ENDIF
 
       ! ----------------------- !
@@ -100,33 +117,40 @@
       zdiff_hc = 0.d0
       zdiff_sc = 0.d0
+
       ! volume variation (calculated with ssh)
-      zdiff_v1 = SUM( surf(:,:) * tmask(:,:,1) * ( sshn(:,:) - ssh_ini(:,:) ) )
+      zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
+
+      ! heat & salt content variation (associated with ssh)
+      IF( .NOT. lk_vvl ) THEN
+         z_ssh_hc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) )
+         z_ssh_sc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) )
+      ENDIF
+
       DO jk = 1, jpkm1
-         ! volume variation (calculated with scale factors)
-         zdiff_v2 = zdiff_v2 + SUM( surf(:,:) * tmask(:,:,jk)   &
+        ! volume variation (calculated with scale factors)
+         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk)   &
             &                       * ( fse3t_n(:,:,jk)         &
             &                           - e3t_ini(:,:,jk) ) )
          ! heat content variation
-         zdiff_hc = zdiff_hc + SUM( surf(:,:) * tmask(:,:,jk)          &
+         zdiff_hc = zdiff_hc + glob_sum( surf(:,:) * tmask(:,:,jk)          &
             &                       * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   &
             &                           - hc_loc_ini(:,:,jk) ) )
          ! salt content variation
-         zdiff_sc = zdiff_sc + SUM( surf(:,:) * tmask(:,:,jk)          &
+         zdiff_sc = zdiff_sc + glob_sum( surf(:,:) * tmask(:,:,jk)          &
             &                       * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   &
             &                           - sc_loc_ini(:,:,jk) ) )
       ENDDO
 
-      IF( lk_mpp ) THEN
-         CALL mpp_sum( zdiff_hc )
-         CALL mpp_sum( zdiff_sc )
-         CALL mpp_sum( zdiff_v1 )
-         CALL mpp_sum( zdiff_v2 )
-      ENDIF
-
       ! Substract forcing from heat content, salt content and volume variations
       zdiff_v1 = zdiff_v1 - frc_v
-      zdiff_v2 = zdiff_v2 - frc_v
+      IF( lk_vvl )   zdiff_v2 = zdiff_v2 - frc_v
       zdiff_hc = zdiff_hc - frc_t
       zdiff_sc = zdiff_sc - frc_s
+      IF( .NOT. lk_vvl ) THEN
+         zdiff_hc1 = zdiff_hc + z_ssh_hc 
+         zdiff_sc1 = zdiff_sc + z_ssh_sc
+         zerr_hc1  = z_ssh_hc - frc_wn_t
+         zerr_sc1  = z_ssh_sc - frc_wn_s
+      ENDIF
       
       ! ----------------------- !
@@ -134,8 +158,15 @@
       ! ----------------------- !
       zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt )
-      WRITE(numhsb , 9020) kt , zdiff_hc / vol_tot , zdiff_hc * fact1  * zdeltat,                                &
-         &                      zdiff_sc / vol_tot , zdiff_sc * fact21 * zdeltat, zdiff_sc * fact22 * zdeltat,   &
-         &                      zdiff_v1           , zdiff_v1 * fact31 * zdeltat, zdiff_v1 * fact32 * zdeltat,   &
-         &                      zdiff_v2           , zdiff_v2 * fact31 * zdeltat, zdiff_v2 * fact32 * zdeltat
+      IF( lk_vvl ) THEN
+         WRITE(numhsb , 9020) kt , zdiff_hc / vol_tot , zdiff_hc * fact1  * zdeltat,                                &
+            &                      zdiff_sc / vol_tot , zdiff_sc * fact21 * zdeltat, zdiff_sc * fact22 * zdeltat,   &
+            &                      zdiff_v1           , zdiff_v1 * fact31 * zdeltat, zdiff_v1 * fact32 * zdeltat,   &
+            &                      zdiff_v2           , zdiff_v2 * fact31 * zdeltat, zdiff_v2 * fact32 * zdeltat
+      ELSE
+         WRITE(numhsb , 9030) kt , zdiff_hc1 / vol_tot , zdiff_hc1 * fact1  * zdeltat,                                &
+            &                      zdiff_sc1 / vol_tot , zdiff_sc1 * fact21 * zdeltat, zdiff_sc1 * fact22 * zdeltat,   &
+            &                      zdiff_v1            , zdiff_v1  * fact31 * zdeltat, zdiff_v1  * fact32 * zdeltat,   &
+            &                      zerr_hc1 / vol_tot  , zerr_sc1 / vol_tot
+      ENDIF
 
       IF ( kt == nitend ) CLOSE( numhsb )
@@ -144,4 +175,5 @@
 
 9020  FORMAT(I5,11D15.7)
+9030  FORMAT(I5,10D15.7)
       !
    END SUBROUTINE dia_hsb
@@ -179,27 +211,18 @@
 
       IF( .NOT. ln_diahsb )   RETURN
+      IF( .NOT. lk_mpp_rep ) &
+        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', &
+             &         ' whereas the global sum to be precise must be done in double precision ',&
+             &         ' please add key_mpp_rep')
 
       ! ------------------- !
       ! 1 - Allocate memory !
       ! ------------------- !
-      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror )
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &
+         &      ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj), &
+         &      e3t_ini(jpi,jpj,jpk)                            , &
+         &      surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror )
       IF( ierror > 0 ) THEN
          CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
-      ENDIF
-      ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' )   ;   RETURN
-      ENDIF
-      ALLOCATE( e3t_ini(jpi,jpj,jpk)   , STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' )      ;   RETURN
-      ENDIF
-      ALLOCATE( surf(jpi,jpj)          , STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate surf' )         ;   RETURN
-      ENDIF
-      ALLOCATE( ssh_ini(jpi,jpj)       , STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' )      ;   RETURN
       ENDIF
 
@@ -214,26 +237,33 @@
       cl_name    = 'heat_salt_volume_budgets.txt'                         ! name of output file
       surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
-      surf_tot  = SUM( surf(:,:) )                                       ! total ocean surface area
+      surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area
       vol_tot   = 0.d0                                                   ! total ocean volume
       DO jk = 1, jpkm1
-         vol_tot  = vol_tot + SUM( surf(:,:) * tmask(:,:,jk)     &
-            &                      * fse3t_n(:,:,jk)         )
+         vol_tot  = vol_tot + glob_sum( surf(:,:) * tmask(:,:,jk)     &
+            &                         * fse3t_n(:,:,jk)         )
       END DO
-      IF( lk_mpp ) THEN 
-         CALL mpp_sum( vol_tot )
-         CALL mpp_sum( surf_tot )
-      ENDIF
 
       CALL ctl_opn( numhsb , cl_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 )
-      !                   12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80
-      WRITE( numhsb, 9010 ) "kt   |     heat content budget     |            salt content budget             ",   &
-         !                                                   123456789012345678901234567890123456789012345 -> 45
-         &                                                  "|            volume budget (ssh)             ",   &
-         !                                                   678901234567890123456789012345678901234567890 -> 45
-         &                                                  "|            volume budget (e3t)             "
-      WRITE( numhsb, 9010 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   &
-         &                                                  "|     [m3]         [mmm/s]          [SV]     ",   &
-         &                                                  "|     [m3]         [mmm/s]          [SV]     "
-
+      IF( lk_vvl ) THEN
+         !                   12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80
+         WRITE( numhsb, 9010 ) "kt   |     heat content budget     |            salt content budget             ",   &
+            !                                                   123456789012345678901234567890123456789012345 -> 45
+            &                                                  "|            volume budget (ssh)             ",   &
+            !                                                   678901234567890123456789012345678901234567890 -> 45
+            &                                                  "|            volume budget (e3t)             "
+         WRITE( numhsb, 9010 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   &
+            &                                                  "|     [m3]         [mmm/s]          [SV]     ",   &
+            &                                                  "|     [m3]         [mmm/s]          [SV]     "
+      ELSE
+         !                   12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80
+         WRITE( numhsb, 9011 ) "kt   |     heat content budget     |            salt content budget             ",   &
+            !                                                   123456789012345678901234567890123456789012345 -> 45
+            &                                                  "|            volume budget (ssh)             ",   &
+            !                                                   678901234567890123456789012345678901234567890 -> 45
+            &                                                  "|  Non conservation due to free surface      "
+         WRITE( numhsb, 9011 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   &
+            &                                                  "|     [m3]         [mmm/s]          [SV]     ",   &
+            &                                                  "|  [heat - C]     [salt - psu]                "
+      ENDIF
       ! --------------- !
       ! 3 - Conversions ! (factors will be multiplied by duration afterwards)
@@ -261,6 +291,13 @@
       frc_t = 0.d0                                           ! heat content   -    -   -    -   
       frc_s = 0.d0                                           ! salt content   -    -   -    -         
+      IF( .NOT. lk_vvl ) THEN
+         ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * ssh_ini(:,:)   ! initial heat content associated with ssh
+         ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * ssh_ini(:,:)   ! initial salt content associated with ssh
+         frc_wn_t = 0.d0
+         frc_wn_s = 0.d0
+      ENDIF
       !
 9010  FORMAT(A80,A45,A45)
+9011  FORMAT(A80,A45,A45)
       !
    END SUBROUTINE dia_hsb_init
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 4177)
@@ -108,15 +108,18 @@
             ncsi1(2)   =  97  ;  ncsj1(2)   = 107
             ncsi2(2)   = 103  ;  ncsj2(2)   = 111
-            ncsir(2,1) = 110  ;  ncsjr(2,1) = 111
-            !                                            ! Black Sea 1 : west part of the Black Sea 
-            ncsnr(3)   = 1    ; ncstt(3)   =   2            !            (ie west of the cyclic b.c.)
-            ncsi1(3)   = 174  ; ncsj1(3)   = 107            ! put in Med Sea
-            ncsi2(3)   = 181  ; ncsj2(3)   = 112
-            ncsir(3,1) = 171  ; ncsjr(3,1) = 106 
-            !                                            ! Black Sea 2 : est part of the Black Sea 
-            ncsnr(4)   =   1  ;  ncstt(4)   =   2           !               (ie est of the cyclic b.c.)
-            ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! put in Med Sea
-            ncsi2(4)   =   6  ;  ncsj2(4)   = 112
-            ncsir(4,1) = 171  ;  ncsjr(4,1) = 106 
+            ncsir(2,1) = 110  ;  ncsjr(2,1) = 111           
+            !                                            ! Black Sea (crossed by the cyclic boundary condition)
+            ncsnr(3:4) =   4  ;  ncstt(3:4) =   2           ! put in Med Sea (north of Aegean Sea)
+            ncsir(3:4,1) = 171;  ncsjr(3:4,1) = 106         !
+            ncsir(3:4,2) = 170;  ncsjr(3:4,2) = 106 
+            ncsir(3:4,3) = 171;  ncsjr(3:4,3) = 105 
+            ncsir(3:4,4) = 170;  ncsjr(3:4,4) = 105 
+            ncsi1(3)   = 174  ;  ncsj1(3)   = 107           ! 1 : west part of the Black Sea      
+            ncsi2(3)   = 181  ;  ncsj2(3)   = 112           !            (ie west of the cyclic b.c.)
+            ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! 2 : east part of the Black Sea 
+            ncsi2(4)   =   6  ;  ncsj2(4)   = 112           !           (ie east of the cyclic b.c.)
+             
+          
+
             !                                        ! =======================
          CASE ( 4 )                                  !  ORCA_R4 configuration
@@ -372,6 +375,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
       !
-      INTEGER  ::   jc, jn      ! dummy loop indices
-      INTEGER  ::   ii, ij      ! temporary integer
+      INTEGER  ::   jc, jn, ji, jj      ! dummy loop indices
       !!----------------------------------------------------------------------
       !
@@ -379,7 +381,9 @@
          IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows
              DO jn = 1, 4
-               ii = mi0( ncsir(jc,jn) )
-               ij = mj0( ncsjr(jc,jn) )
-               p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0_wp )
+                DO jj =    mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) )
+                   DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) )
+                      p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp )
+                   END DO
+                END DO
             END DO 
          ENDIF 
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90	(revision 4177)
@@ -238,4 +238,9 @@
                nday_year = 1
                nsec_year = ndt05
+               IF( nsec1jan000 >= 2 * (2**30 - nsecd * nyear_len(1) / 2 ) ) THEN   ! test integer 4 max value
+                  CALL ctl_stop( 'The number of seconds between Jan. 1st 00h of nit000 year and Jan. 1st 00h ',   &
+                     &           'of the current year is exceeding the INTEGER 4 max VALUE: 2^31-1 -> 68.09 years in seconds', &
+                     & 'You must do a restart at higher frequency (or remove this STOP and recompile everything in I8)' )
+               ENDIF
                nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1)
                IF( nleapy == 1 )   CALL day_mth
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 4177)
@@ -1102,9 +1102,6 @@
       INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
       REAL(wp) ::   zrmax, ztaper   ! temporary scalars
-      REAL(wp) ::   zrfact   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
-
-      !
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, zri, zrj, zhbat
+      !
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
 
       NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
@@ -1114,7 +1111,6 @@
       IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
       !
-      CALL wrk_alloc( jpi, jpj,      ztmpi1, ztmpi2, ztmpj1, ztmpj2         )
-      CALL wrk_alloc( jpi, jpj,      zenv, zri, zrj, zhbat     )
-     !
+      CALL wrk_alloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
+      !
       REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters
       READ  ( numnam, namzgr_sco )
@@ -1163,11 +1159,9 @@
       !                                        ! =============================
       ! use r-value to create hybrid coordinates
-!     DO jj = 1, jpj
-!        DO ji = 1, jpi
-!           zenv(ji,jj) = MAX( bathy(ji,jj), 0._wp )
-!        END DO
-!     END DO
-!     CALL lbc_lnk( zenv, 'T', 1._wp )
-      zenv(:,:) = bathy(:,:)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
+         END DO
+      END DO
       ! 
       ! Smooth the bathymetry (if required)
@@ -1177,78 +1171,87 @@
       jl = 0
       zrmax = 1._wp
-      !     
-      ! set scaling factor used in reducing vertical gradients
-      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 
-      !
-      ! initialise temporary evelope depth arrays
-      ztmpi1(:,:) = zenv(:,:)
-      ztmpi2(:,:) = zenv(:,:)
-      ztmpj1(:,:) = zenv(:,:)
-      ztmpj2(:,:) = zenv(:,:)
-      !
-      ! initialise temporary r-value arrays
-      zri(:,:) = 1._wp
-      zrj(:,:) = 1._wp
-      !                                                            ! ================ !
-      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
-         !                                                         ! ================ !
+      !                                                     ! ================ !
+      DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )         !  Iterative loop  !
+         !                                                  ! ================ !
          jl = jl + 1
          zrmax = 0._wp
-         ! we set zrmax from previous r-values (zri abd zrj) first
-         ! if set after current r-value calculation (as previously)
-         ! we could exit DO WHILE prematurely before checking r-value
-         ! of current zenv
-         DO jj = 1, nlcj
-            DO ji = 1, nlci
-               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
-            END DO
-         END DO
-         zri(:,:) = 0._wp
-         zrj(:,:) = 0._wp
+         zmsk(:,:) = 0._wp
          DO jj = 1, nlcj
             DO ji = 1, nlci
                iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
                ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
-               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
-                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
-               END IF
-               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
-                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
-               END IF
-               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
-               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact 
-               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
-               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
+               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
+               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
+               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
+               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
+               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1._wp
+               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
+               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp
             END DO
          END DO
          IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
+         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
+         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1._wp )
+         DO jj = 1, nlcj
+            DO ji = 1, nlci
+                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
+            END DO
+         END DO
          !
-         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
+         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
          !
          DO jj = 1, nlcj
             DO ji = 1, nlci
-               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
+               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
+               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
+               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
+               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
+               IF( zmsk(ji,jj) == 1._wp ) THEN
+                  ztmp(ji,jj) =   (                                                                                   &
+             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
+             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
+             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
+             &                    ) / (                                                                               &
+             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
+             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   &
+             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
+             &                        )
+               ENDIF
             END DO
          END DO
          !
-         CALL lbc_lnk( zenv, 'T', 1._wp )
+         DO jj = 1, nlcj
+            DO ji = 1, nlci
+               IF( zmsk(ji,jj) == 1._wp )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
+            END DO
+         END DO
+         !
+         ! Apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
+         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1._wp )
+         DO jj = 1, nlcj
+            DO ji = 1, nlci
+               IF( zenv(ji,jj) == 0._wp )   zenv(ji,jj) = ztmp(ji,jj)
+            END DO
+         END DO
          !                                                  ! ================ !
       END DO                                                !     End loop     !
       !                                                     ! ================ !
       !
-!     DO jj = 1, jpj
-!        DO ji = 1, jpi
-!           zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale values
-!        END DO
-!     END DO
+      ! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions
+      DO ji = nlci+1, jpi 
+         zenv(ji,1:nlcj) = zenv(nlci,1:nlcj)
+      END DO
+      !
+      DO jj = nlcj+1, jpj
+         zenv(:,jj) = zenv(:,nlcj)
+      END DO
       !
       ! Envelope bathymetry saved in hbatt
       hbatt(:,:) = zenv(:,:) 
-
       IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
          CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )
+               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
                hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
             END DO
@@ -1365,11 +1368,11 @@
       fsde3w(:,:,:) = gdep3w(:,:,:)
       !
-      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0
-      where (e3u   (:,:,:).eq.0.0)  e3u(:,:,:) = 1.0
-      where (e3v   (:,:,:).eq.0.0)  e3v(:,:,:) = 1.0
-      where (e3f   (:,:,:).eq.0.0)  e3f(:,:,:) = 1.0
-      where (e3w   (:,:,:).eq.0.0)  e3w(:,:,:) = 1.0
-      where (e3uw  (:,:,:).eq.0.0)  e3uw(:,:,:) = 1.0
-      where (e3vw  (:,:,:).eq.0.0)  e3vw(:,:,:) = 1.0
+      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1._wp
+      where (e3u   (:,:,:).eq.0.0)  e3u(:,:,:) = 1._wp
+      where (e3v   (:,:,:).eq.0.0)  e3v(:,:,:) = 1._wp
+      where (e3f   (:,:,:).eq.0.0)  e3f(:,:,:) = 1._wp
+      where (e3w   (:,:,:).eq.0.0)  e3w(:,:,:) = 1._wp
+      where (e3uw  (:,:,:).eq.0.0)  e3uw(:,:,:) = 1._wp
+      where (e3vw  (:,:,:).eq.0.0)  e3vw(:,:,:) = 1._wp
 
 #if defined key_agrif
@@ -1519,5 +1522,6 @@
       END DO
       !
-      CALL wrk_dealloc( jpi, jpj,      zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat                           )      !
+      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
+      !
       IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
       !
@@ -1748,9 +1752,4 @@
       ENDDO
       !
-      CALL lbc_lnk(e3t ,'T',1.) ; CALL lbc_lnk(e3u ,'T',1.)
-      CALL lbc_lnk(e3v ,'T',1.) ; CALL lbc_lnk(e3f ,'T',1.)
-      CALL lbc_lnk(e3w ,'T',1.)
-      CALL lbc_lnk(e3uw,'T',1.) ; CALL lbc_lnk(e3vw,'T',1.)
-      !
       !                                               ! =============
 
@@ -1849,5 +1848,5 @@
       !!----------------------------------------------------------------------
       !
-      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
+      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1,wp) + rn_thetb )  )   &
          &     - TANH( rn_thetb * rn_theta                                )  )   &
          & * (   COSH( rn_theta                           )                      &
@@ -1875,8 +1874,8 @@
       !
       IF ( rn_theta == 0 ) then      ! uniform sigma
-         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
+         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1,wp )
       ELSE                        ! stretched sigma
-         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
-            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
+         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1,wp)) ) ) / SINH( rn_theta )              &
+            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1,wp)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
             &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
       ENDIF
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 4177)
@@ -109,10 +109,8 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zub, zvb
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt')
       !
-      CALL wrk_alloc( jpi,jpj,jpk, zub, zvb )
       !
       IF( kt == nit000 ) THEN
@@ -213,6 +211,6 @@
          DO jk = 1, jpkm1
             DO ji = 1, jpij
-               spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk)
-               spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk)
+               spgu(ji,1) = spgu(ji,1) + fse3u_a(ji,1,jk) * ua(ji,1,jk)
+               spgv(ji,1) = spgv(ji,1) + fse3v_a(ji,1,jk) * va(ji,1,jk)
             END DO
          END DO
@@ -221,6 +219,6 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)
-                  spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)
+                  spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)
+                  spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)
                END DO
             END DO
@@ -360,5 +358,4 @@
       IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
       !
-      CALL wrk_dealloc( jpi,jpj,jpk, zub, zvb )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt')
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 4177)
@@ -2181,9 +2181,9 @@
 !!gm Remark : this is very time consumming!!!
       !                                         ! ------------------------ !
-            IF( ijpt0 > ijpt1 .OR. iipt0 > iipt1 ) THEN
+        IF(((nbondi .ne. 0) .AND. (ktype .eq. 2)) .OR. ((nbondj .ne. 0) .AND. (ktype .eq. 1))) THEN
             ! there is nothing to be migrated
-               lmigr = .FALSE.
+              lmigr = .TRUE.
             ELSE
-              lmigr = .TRUE.
+              lmigr = .FALSE.
             ENDIF
 
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 4177)
@@ -187,5 +187,5 @@
          &      gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr )
       IF(lk_mpp)   CALL mpp_sum( ierr )
-      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' )
+      IF( ierr /= 0 )   CALL ctl_stop('angle: unable to allocate arrays' )
 
       ! ============================= !
@@ -361,5 +361,5 @@
             &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
          IF( lk_mpp    )   CALL mpp_sum( ierr )
-         IF( ierr /= 0 )   CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' )
+         IF( ierr /= 0 )   CALL ctl_stop('geo2oce: unable to allocate arrays' )
       ENDIF
 
@@ -438,9 +438,9 @@
       !!----------------------------------------------------------------------
 
-      IF( ALLOCATED( gsinlon ) ) THEN
+      IF( .NOT. ALLOCATED( gsinlon ) ) THEN
          ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   &
             &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
          IF( lk_mpp    )   CALL mpp_sum( ierr )
-         IF( ierr /= 0 )   CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' )
+         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' )
       ENDIF
 
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 4177)
@@ -388,5 +388,5 @@
       !
       IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
-         srcv(jpr_itz1:jpr_itz2)%laction = .FALSE.    ! ice components not received (itx1 and ity1 used later)
+         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
          srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
          srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
@@ -407,5 +407,7 @@
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
       CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
-      CASE( 'conservative'  )   ;   srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
+      CASE( 'conservative'  )
+         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
+         IF ( k_ice <= 1 )  srcv(jpr_ivep)%laction = .FALSE.
       CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
       CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
@@ -465,13 +467,4 @@
          CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
       !                                                      ! ------------------------- !
-      !                                                      !    Ice Qsr penetration    !   
-      !                                                      ! ------------------------- !
-      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer 
-      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
-      ! Coupled case: since cloud cover is not received from atmosphere 
-      !               ===> defined as constant value -> definition done in sbc_cpl_init
-      fr1_i0(:,:) = 0.18
-      fr2_i0(:,:) = 0.82
-      !                                                      ! ------------------------- !
       !                                                      !      10m wind module      !   
       !                                                      ! ------------------------- !
@@ -508,4 +501,9 @@
       ! Allocate taum part of frcv which is used even when not received as coupling field
       IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jn)%nct) )
+      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
+      IF( k_ice /= 0 ) THEN
+         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jn)%nct) )
+         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jn)%nct) )
+      END IF
 
       ! ================================ !
@@ -1329,4 +1327,13 @@
       END SELECT
 
+      !    Ice Qsr penetration used (only?)in lim2 or lim3 
+      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer 
+      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
+      ! Coupled case: since cloud cover is not received from atmosphere 
+      !               ===> defined as constant value -> definition done in sbc_cpl_init
+      fr1_i0(:,:) = 0.18
+      fr2_i0(:,:) = 0.82
+
+
       CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
       !
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4177)
@@ -221,4 +221,6 @@
       ENDIF
       !
+                          CALL sbc_ssm_init               ! Sea-surface mean fields initialisation
+      !
       IF( ln_ssr      )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation
       !
Index: branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 4176)
+++ branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 4177)
@@ -675,5 +675,5 @@
 
 
-   FUNCTION tfreez( psal ) RESULT( ptf )
+   FUNCTION tfreez( psal, pdep ) RESULT( ptf )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE eos_init  ***
@@ -688,4 +688,5 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [decibars]
       ! Leave result array automatic rather than making explicitly allocated
       REAL(wp), DIMENSION(jpi,jpj)                ::   ptf    ! freezing temperature [Celcius]
@@ -694,4 +695,7 @@
       ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
          &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
+      IF ( PRESENT( pdep ) ) THEN   
+         ptf(:,:) = ptf(:,:) - 7.53e-4_wp * pdep(:,:)
+      ENDIF
       !
    END FUNCTION tfreez
