MODULE bulk
   !!======================================================================
   !!                           ***  bulk  ***
   !!======================================================================
#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
   !!----------------------------------------------------------------------
   !!   'key_flx_bulk_monthly'                        monthly bulk formulea
   !!   'key_flx_bulk_daily'                          monthly bulk formulea
   !!----------------------------------------------------------------------
   !!   bulk          : computation of fluxes using bulk formulation
   !!----------------------------------------------------------------------
   !! * Modules used   
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE ice_oce         ! bulk variable  
   USE ocfzpt          ! ocean freezing point
   USE flxblk          ! bulk formulae
   USE blk_oce         ! bulk variable 
   USE flx_oce
   USE taumod
   USE phycst          ! physical constants
   USE in_out_manager  ! I/O manager

   IMPLICIT NONE
   PRIVATE

   !! * Routine accessibility
   PUBLIC blk        ! called by flx.F90   
   !!----------------------------------------------------------------------
   !!   OPA 9.0 , LODYC-IPSL  (2003)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE blk( kt )
      !!---------------------------------------------------------------------
      !!                    ***  ROUTINE blk  ***
      !!        
      !! ** Purpose :   provide the heat fluxes on ice and ocean 
      !!                using bulk formulation
      !!
      !! History :
      !!   9.0  !  03-11  (C. Ethe and G. Madec)  F90: Free form and MODULE 
      !!----------------------------------------------------------------------
      !! * arguments
      INTEGER, INTENT( in  ) ::   kt   ! ocean time step

      !! * Local declarations    
      REAL(wp), DIMENSION(jpi,jpj) ::   &
         zsst , zsss
# if ! defined key_ice_lim
      INTEGER  ::   ji, jj         ! dummy loop indices  
      REAL(wp) ::   ztgel, zicopa
# endif
      !!---------------------------------------------------------------------

     ! Initialisation
     IF( kt == nit000) THEN
      ! computation of rdtbs2
        IF( nacc == 1 ) THEN
           rdtbs2 = nfbulk * rdtmin * 0.5
        ELSE
           rdtbs2 = nfbulk * rdt * 0.5
        ENDIF
        IF ( .NOT.ln_rstart ) THEN
           gsst(:,:) =  ( nfbulk - 1 ) * ( tn(:,:,1) + rt0 )
           gsss(:,:) =  ( nfbulk - 1 ) *   sn(:,:,1)
        ENDIF
     ENDIF

# if ! defined key_ice_lim
      ! opa model ice freeze()      
      DO jj = 1, jpj
         DO ji = 1, jpi
            ztgel  = fzptn(ji,jj)
            zicopa = tmask(ji,jj,1)
            IF( tn(ji,jj,1) >= ztgel ) zicopa = 0.
            freeze(ji,jj) = zicopa
         END DO
      END DO
# endif

      gsst(:,:) = gsst(:,:) + tn(:,:,1) + rt0  
      gsss(:,:) = gsss(:,:) + sn(:,:,1)

      !  Computation of the fluxes       
      IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN

!i       zsst(:,:) = gsst(:,:) / REAL( nfbulk ) 
         zsst(:,:) = gsst(:,:) / REAL( nfbulk ) * tmask(:,:,1)
         zsss(:,:) = gsss(:,:) / REAL( nfbulk )
         CALL flx_blk( kt, zsst )    
     
         gsst(:,:) = 0.    
         gsss(:,:) = 0.  

# if ! defined key_ice_lim
         IF(l_ctl) THEN         ! print mean trends (used for debugging)
            WRITE(numout,*) ' Forcings '
            WRITE(numout,*) ' qsr_oce  : ', SUM( qsr_oce (:,:) * tmask(:,:,1) )
            WRITE(numout,*) ' qsr_ice  : ', SUM( qsr_ice (:,:) * tmask(:,:,1) )
            WRITE(numout,*) ' qnsr_oce : ', SUM( qnsr_oce(:,:) * tmask(:,:,1) )
            WRITE(numout,*) ' qnsr_ice : ', SUM( qnsr_ice(:,:) * tmask(:,:,1) )
            WRITE(numout,*) ' evap     : ', SUM( evap    (:,:) * tmask(:,:,1) )
            WRITE(numout,*) ' precip   : ', SUM( tprecip (:,:) * tmask(:,:,1) ) / rday
            WRITE(numout,*) ' Snow     : ', SUM( sprecip (:,:) * tmask(:,:,1) ) / rday
            WRITE(numout,*) ' u-stress : ', SUM( taux    (:,:) * umask(:,:,1) )
            WRITE(numout,*) ' v-stress : ', SUM( tauy    (:,:) * vmask(:,:,1) )
            WRITE(numout,*) ' sst      : ', SUM( zsst    (:,:) * tmask(:,:,1) )
            WRITE(numout,*) ' sss      : ', SUM( zsss    (:,:) * tmask(:,:,1) )
            WRITE(numout,*)
         ENDIF
# endif   
      ENDIF
 
   END SUBROUTINE blk

#else
   !!----------------------------------------------------------------------
   !!   Dummy module :                                     NO bulk formulea
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE blk( kt )          ! Dummy routine
      WRITE(*,*) 'blk: You should not see this print! error? ', kt
   END SUBROUTINE blk
#endif
 
   !!======================================================================
END MODULE bulk