!!----------------------------------------------------------------------
   !!                       ***  ocesbc_cpl_ice.h90  ***
   !!----------------------------------------------------------------------

   SUBROUTINE oce_sbc(kt)
      !!---------------------------------------------------------------------
      !!                 ***  ROUTINE oce_sbc  ***
      !!                    
      !! ** Purpose :   Ocean surface boundaries conditions with 
      !!        Louvain la Neuve Sea Ice Model in coupled mode
      !!
      !! History :
      !!   1.0  !  00-10  (O. Marti)  Original code
      !!   2.0  !  02-12  (G. Madec)  F90: Free form and module
      !!----------------------------------------------------------------------
      !! * Local declarations
      INTEGER ::   ji, jj                   ! dummy loop indices
      REAL(wp) ::   ztx, ztaux, zty, ztauy, ztrp, zsrp
      REAL(wp) ::   ztdta, ztgel, zqrp

      !!----------------------------------------------------------------------
      !!  OPA 8.5, LODYC-IPSL (2002)
      !!----------------------------------------------------------------------
 
      !!----------------------------------------------------------------------        
      ! 1. initialization to zero at kt = nit000
      ! ---------------------------------------
      
      IF( kt == nit000 ) THEN     
         qsr   (:,:) = 0.e0
         freeze(:,:) = 0.e0
         qt    (:,:) = 0.e0
         q     (:,:) = 0.e0
         qrp   (:,:) = 0.e0
         emp   (:,:) = 0.e0
         emps  (:,:) = 0.e0
         erp   (:,:) = 0.e0
#if defined key_dynspg_fsc 
         dmp   (:,:) = 0.e0
#endif
      ENDIF

      IF( MOD( kt-1, nfice ) == 0 ) THEN 

         CALL oce_sbc_dmp   ! Computation of internal and evaporation damping terms       

         ! Surface heat flux (W/m2)
         ! -----------------------

         ! restoring heat flux
         DO jj = 1, jpj
            DO ji = 1, jpi
               ztgel = fzptn(ji,jj)
#if defined key_dtasst
               ztdta = MAX( sst(ji,jj),    ztgel )
#else
               ztdta = MAX( t_dta(ji,jj,1), ztgel )
#endif
               zqrp = ztrp * ( tb(ji,jj,1) - ztdta )

               qrp(ji,jj) = (1.0-freeze(ji,jj) ) * zqrp
            END DO
         END DO

         ! non solar heat flux + solar flux + restoring
         q  (:,:) = fnsolar(:,:) + fsolar(:,:) + qrp(:,:)
         qt (:,:) = q(:,:)

         ! solar flux
         qsr(:,:) = fsolar(:,:)

#if defined key_dynspg_fsc  
         ! total concentration/dilution effect (use on SSS)
         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)

         ! total volume flux (use on sea-surface height)
         emp (:,:) = fmass(:,:)  -  dmp(:,:) + runoff(:,:) + erp(:,:)

#else
         ! Rigid-lid (emp=emps=E-P-R+Erp) 
         ! freshwater flux
         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
         emp (:,:) = emps(:,:)

#endif

         DO jj = 1, jpjm1
            DO ji = 1, fs_jpim1   ! vertor opt.
               ztx   = 0.5 * ( freeze(ji+1,jj) + freeze(ji+1,jj+1) )
               ztaux = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) )
               taux(ji,jj) = (1.0-ztx) * taux(ji,jj) + ztx * ztaux

               zty   = 0.5 * ( freeze(ji,jj+1) + freeze(ji+1,jj+1) )
               ztauy = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) )
               tauy(ji,jj) = (1.0-zty) * tauy(ji,jj) + zty * ztauy
            END DO
         END DO
         CALL lbc_lnk( taux, 'U', -1. )
         CALL lbc_lnk( tauy, 'V', -1. )    

         ! Re-initialization of fluxes
         sst_io(:,:) = 0.0
         sss_io(:,:) = 0.0
         u_io  (:,:) = 0.0
         v_io  (:,:) = 0.0
         gtaux (:,:) = 0.
         gtauy (:,:) = 0.

      ENDIF

   END SUBROUTINE oce_sbc