Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 3413)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 3414)
@@ -158,5 +158,4 @@
          &                       rdm_ice(jpi,jpj) , rdq_ice(jpi,jpj) ,     &
          &     qldif (jpi,jpj) , qcmif  (jpi,jpj) , fdtcn  (jpi,jpj) ,     &
-         &     qldif (jpi,jpj) , qcmif  (jpi,jpj) , fdtcn  (jpi,jpj) ,     &
          &     qdtcn (jpi,jpj) , thcm   (jpi,jpj)                    , STAT=ierr(4) )
 
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 3413)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 3414)
@@ -30,4 +30,5 @@
    USE in_out_manager ! I/O manager
    USE prtctl         ! Print control
+   USE oce, ONLY : snwice_mass, snwice_mass_b
 
    IMPLICIT NONE
@@ -80,8 +81,9 @@
       REAL(wp) ::   zs21_11, zs21_12, zs21_21, zs21_22
       REAL(wp) ::   zs22_11, zs22_12, zs22_21, zs22_22
+      REAL(wp) ::   zintb, zintn
       REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld, zmass, zcorl
       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct, za2ct, zresr
       REAL(wp), POINTER, DIMENSION(:,:) ::   zc1u, zc1v, zc2u, zc2v
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zsang
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zsang, zpice
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu0, zv0
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_n, zv_n
@@ -93,5 +95,5 @@
       
       CALL wrk_alloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr )
-      CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang )
+      CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice )
       CALL wrk_alloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 )
       CALL wrk_alloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 )
@@ -129,4 +131,20 @@
 !i    zviszeta(:,jpj+1) = 0._wp    ;    zviseta(:,jpj+1) = 0._wp
 
+      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
+          !
+          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
+         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
+          !
+          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
+         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
+          !
+         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
+          !
+         !
+      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
+         zpice(:,:) = ssh_m(:,:)
+      ENDIF
 
       ! Ice mass, ice strength, and wind stress at the center            |
@@ -196,8 +214,8 @@
 
             ! Gradient of the sea surface height
-            zgsshx =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji-1,jj  ))/e1u(ji-1,jj  )   &
-               &       +  (ssh_m(ji  ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp
-            zgsshy =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji  ,jj-1))/e2v(ji  ,jj-1)   &
-               &       +  (ssh_m(ji-1,jj  ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp
+            zgsshx =  (   (zpice(ji  ,jj  ) - zpice(ji-1,jj  ))/e1u(ji-1,jj  )   &
+               &       +  (zpice(ji  ,jj-1) - zpice(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp
+            zgsshy =  (   (zpice(ji  ,jj  ) - zpice(ji  ,jj-1))/e2v(ji  ,jj-1)   &
+               &       +  (zpice(ji-1,jj  ) - zpice(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp
 
             ! Computation of the velocity field taking into account the ice-ice interaction.                                 
@@ -575,5 +593,5 @@
 
       CALL wrk_dealloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr )
-      CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang )
+      CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice )
       CALL wrk_dealloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 )
       CALL wrk_dealloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 )
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3413)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3414)
@@ -470,5 +470,5 @@
       ENDIF
       IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart : 
-         &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh belew sea-ice area
+         &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area
          sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
          sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 3413)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 3414)
@@ -35,4 +35,5 @@
    USE dom_ice_2        ! LIM2: ice domain
 #endif
+   USE oce, ONLY :  snwice_mass, snwice_mass_b
 
    IMPLICIT NONE
@@ -124,4 +125,5 @@
       REAL(wp) ::   zindb         ! ice (1) or not (0)      
       REAL(wp) ::   zdummy        ! dummy argument
+      REAL(wp) ::   zintb, zintn  ! dummy argument
 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength
@@ -144,4 +146,7 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope:
+                                                              !   ocean surface (ssh_m) if ice is not embedded
+                                                              !   ice top surface if ice is embedded
       
       !!-------------------------------------------------------------------
@@ -150,5 +155,5 @@
       CALL wrk_alloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
       CALL wrk_alloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds          )
-      CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr                         )
+      CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
 
 #if  defined key_lim2 && ! defined key_lim2_vp
@@ -231,4 +236,20 @@
       !  v_oce2: ocean v component on v points                        
 
+      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
+          !                                            
+          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
+         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
+          !
+          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
+         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
+          !
+         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0
+          !
+      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
+         zpice(:,:) = ssh_m(:,:)
+      ENDIF
+
       DO jj = k_j1+1, k_jpj-1
          DO ji = fs_2, fs_jpim1
@@ -273,6 +294,6 @@
             ! include it later
 
-            zdsshx =  ( ssh_m(ji+1,jj) - ssh_m(ji,jj) ) / e1u(ji,jj)
-            zdsshy =  ( ssh_m(ji,jj+1) - ssh_m(ji,jj) ) / e2v(ji,jj)
+            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
+            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
 
             za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
@@ -746,5 +767,5 @@
       CALL wrk_dealloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                )
       CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds          )
-      CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr                         )
+      CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
 
    END SUBROUTINE lim_rhg
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 3413)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 3414)
@@ -255,5 +255,5 @@
       wndm(:,:)        = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1)    ! wind speed - needed for gas exchange
       emp (:,:)        = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
-      emps(:,:)        = emp(:,:) 
+      sfx (:,:)        = emp(:,:) 
       fr_i(:,:)        = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)     ! Sea-ice fraction
       qsr (:,:)        = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation
@@ -302,5 +302,5 @@
          CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=emps             , clinfo1=' emps    - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=sfx              , clinfo1=' sfx     - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 )
          CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 )
