!!----------------------------------------------------------------------
   !!                     ***  ocesbc_forced_noice.h90  ***
   !!----------------------------------------------------------------------

   SUBROUTINE oce_sbc(kt)
      !!---------------------------------------------------------------------
      !!                   ***    ROUTINE oce_sbc ***
      !!                    
      !! ** Purpose : Ocean surface boundary conditions 
      !!        in forced mode using flux formulation or bulk formulation
      !!       
      !!
      !! 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
      !!----------------------------------------------------------------------
      !! * Modules used
      USE daymod
#if ! defined key_dtasst
      USE dtasst, ONLY : rclice
#endif

#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
      USE blk_oce
#endif

#if defined key_flx_forced_daily
      USE flx_oce
#endif
      !! * arguments
      INTEGER, INTENT( in  ) ::   kt   ! ocean time step

      !! * local declarations
      INTEGER ::   ji, jj        ! dummy loop arguments
      INTEGER ::   i15, ifreq             !      
      REAL(wp) ::  zxy
      REAL(wp) ::  zsice, zqri, zqrp, ztdta, zqrj
      REAL(wp) ::  zq, zqi, zhemis, ztrp
      REAL(wp), DIMENSION(jpi,jpj) :: zeri, zerps, ziclim
      REAL(wp), DIMENSION(jpi,jpj) :: zqt, zqsr, zemp  

      !!----------------------------------------------------------------------
      !!  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 defined key_flx_bulk_monthly || defined key_flx_bulk_daily
      ifreq      = nfbulk
      zqt (:,:)  = qsr_oce(:,:) + qnsr_oce(:,:)
      zqsr(:,:)  = qsr_oce(:,:)
      zemp(:,:)  = evap(:,:) - tprecip(:,:)
#endif 

#if defined key_flx_forced_daily
      ifreq      = 1
      zqt (:,:)  = p_qt (:,:)
      zqsr(:,:)  = p_qsr(:,:)
      zemp(:,:)  = p_emp(:,:)
#endif 

      IF( MOD( kt-1, ifreq) == 0 ) THEN
         ! Computation of internal and evaporation damping terms       
         CALL oce_sbc_dmp

         ztrp = -40.             ! restoring terme for temperature (w/m2/k)   
         zsice = - 0.04 / 0.8    ! ratio of isohaline compressibility over isotherme compressibility 
                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
         ! Flux computation
         DO jj = 1, jpj
            DO ji = 1, jpi      
               ! climatological ice 
#if defined key_dtasst 
               ziclim(ji,jj) = FLOAT( NINT( rclice(ji,jj,1) ) )
#else
               ! tested only with key key_dtasst (A. Lazar 07/2001)
               ! this loop in CASE of interpolation of monthly rclice
               i15           = INT( 2.* FLOAT(nday) / (FLOAT( nobis(nmonth) ) + 0.5) )
               zxy           = FLOAT(nday) / FLOAT(nobis(nmonth)) + 0.5 - i15
               ziclim(ji,jj) = FLOAT( NINT( (1-zxy) * rclice(ji,jj,1) + zxy  * rclice(ji,jj,2) ) )
#endif

               ! avoid surfreezing point            
               tn(ji,jj,1) = MAX( tn(ji,jj,1), fzptn(ji,jj) )

               ! hemisphere indicator (=1 north, =-1 south)           
               zhemis = FLOAT( isign(1, mjg(jj) - (jpjdta/2+1) ) )

               ! restoring temperature (ztdta >= to local freezing temperature)            
#if defined key_dtasst
               ztdta = MAX( sst(ji,jj),    fzptn(ji,jj) )
#else
               ztdta = MAX( t_dta(ji,jj,1), fzptn(ji,jj) )
#endif

               ! a) net downward radiative flux qsr()           
               qsr(ji,jj) = zqsr(ji,jj) * tmask(ji,jj,1)

               ! b) heat flux damping term qrp()
               ! - gamma*(t-tlevitus) if no  climatological ice (ziclim=0)
               ! - gamma*(t-(tgel-1.))  if climatological ice and no opa ice   (ziclim=1 zicopa=0)
               ! - gamma*min(0,t-tgel) if climatological and opa ice (ziclim=1 zicopa=1)

               zqrp = ztrp * ( tb(ji,jj,1) - ztdta )
               zqri = ztrp * ( tb(ji,jj,1) - ( fzptn(ji,jj) - 1.) )
               zqrj = ztrp * MIN( 0., tb(ji,jj,1) - fzptn(ji,jj) )
               qrp(ji,jj) = ( (1. - ziclim(ji,jj)) * zqrp   &
                  + ziclim(ji,jj)  * ( ( 1 - freeze(ji,jj)) * zqri   &
                  + freeze(ji,jj)  * zqrj ) ) * tmask(ji,jj,1)

               ! c) net downward heat flux q() = q0 + qrp()
               ! for q0
               ! - ECMWF fluxes if no climatological ice      (ziclim=0)
               ! - qrp if climatological ice and no opa ice   (ziclim=1 zicopa=0)
               ! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if climatological and opa ice 
               !                                              (ziclim=1 zicopa=1)
               zq  = zqt(ji,jj)
               zqi = -3. + zhemis
               qt (ji,jj) = ( (1.-ziclim(ji,jj)) * zq   &
                  +ziclim(ji,jj)  * freeze(ji,jj) * zqi )   &
                  * tmask(ji,jj,1)   &
                  + qrp(ji,jj)
               q  (ji,jj) = qt (ji,jj)

            END DO
         END DO

#if defined key_dynspg_fsc 
         ! Free-surface

         ! Water flux for zero buoyancy flux if no opa ice and ice clim
         zeri(:,:) = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )

         ! Contribution to sea level:
         ! net upward water flux emp() = e-p + runoff() + erp() + dmp + empold
         emp (:,:) = zemp(:,:)     &   ! e-p data
                   + runoff(:,:)   &   ! runoff data
                   + erp(:,:)      &   ! restoring term to SSS data
                   + dmp(:,:)      &   ! freshwater flux associated with internal damping
                   + empold            ! domain averaged annual mean correction

         ! Contribution to salinity:
         ! net upward water flux emps() = e-p + runoff() + erp() + zerps + empold
         emps(:,:) = zemp(:,:)     &
                   + runoff(:,:)   &
                   + erp(:,:)      &
                   + zerps(:,:)    &
                   + empold 
#else
         ! Rigid-lid (emp=emps=E-P-R+Erp)
         ! freshwater flux
         zeri(:,:)  = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
         zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
         emps (:,:) = zemp(:,:)     &
                    + runoff(:,:)   &
                    + erp(:,:)      &
                    + zerps(:,:)
         emp (:,:) = emps(:,:)
#endif  


         ! Boundary condition on emp for free surface option
         ! -------------------------------------------------
         CALL lbc_lnk( emp, 'T', 1. )
      
      ENDIF

   END SUBROUTINE oce_sbc