!!----------------------------------------------------------------------
   !!                     ***  ocesbc_forced_ice.h90  ***
   !!----------------------------------------------------------------------

   SUBROUTINE oce_sbc (kt)
      !!---------------------------------------------------------------------
      !!                   ***    ROUTINE oce_sbc ***
      !!                    
      !! ** Purpose : Ocean surface boundary conditions with 
      !!        Louvain la Neuve Sea Ice Model in forced mode
      !!        using 'clio' bulk formulea
      !!
      !! History :
      !!   1.0  !  99-11  (M. Imbard)  Original code
      !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
      !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
      !!----------------------------------------------------------------------
      !! * arguments
      INTEGER, INTENT( in  ) ::   kt   ! ocean time step

      !! * Local declarations
      INTEGER ::   ji, jj                   ! dummy loop indices
      REAL ztx, ztaux, zty, ztauy

      !!----------------------------------------------------------------------
      !!  OPA 8.5, LODYC-IPSL (2002)
      !!----------------------------------------------------------------------
 

      ! 1. initialization to zero at kt = nit000
      ! ---------------------------------------
      
      IF( kt == nit000 ) THEN     
         qsr    (:,:) = 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

         ! Computation of internal and evaporation damping terms       
         CALL oce_sbc_dmp

         ! Surface Ocean fluxes
         ! ====================

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

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

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

         ! total volume flux (use on sea-surface height)
         emp (:,:) = fmass(:,:) -   dmp(:,:) + runoff(:,:) + erp(:,:) + empold      
#else
         ! Rigid-lid (emp=emps=E-P-R+Erp)
         emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)     ! freshwater flux
         emp (:,:) = emps(:,:)

#endif

         ! Surface stress
         ! --------------

         ! update the stress beloww sea-ice area
         DO jj = 1, jpjm1
            DO ji = 1, fs_jpim1   ! vertor opt.
               ztx         = MAX( freezn(ji,jj), freezn(ji,jj+1) )   ! ice/ocean indicator at U- and V-points
               zty         = MAX( freezn(ji,jj), freezn(ji+1,jj) )
               ztaux       = 0.5 *( ftaux(ji+1,jj) + ftaux(ji+1,jj+1) ) ! ice-ocean stress at U- and V-points
               ztauy       = 0.5 *( ftauy(ji,jj+1) + ftauy(ji+1,jj+1) )
               taux(ji,jj) = (1.-ztx) * taux(ji,jj) + ztx * ztaux    ! stress at the ocean surface
               tauy(ji,jj) = (1.-zty) * tauy(ji,jj) + zty * ztauy
            END DO
         END DO

         ! boundary condition on the stress (taux,tauy)
         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

      ENDIF

   END SUBROUTINE oce_sbc