SUBROUTINE ice_th(nlay_i,nlay_s)

        !!------------------------------------------------------------------
        !!                ***         ROUTINE ice_th       ***
        !! ** Purpose :
        !!           Ice thermodynamics 
        !! ** Method  :
        !!    This routine calls the thermodynamic and biological routines
        !!
        !! ** Arguments :
        !!           nlay_i, nlay_s
        !!
        !! ** Inputs / Ouputs : (global commons)
        !!
        !! ** External : 
        !!
        !! ** References : Vancoppenolle et al., JGR 2007
        !!
        !! ** History :
        !!       (1) CLIO, Goosse and Fichefet, JGR, 1999.
        !!       (2) LIM-1D, Vancoppenolle et al., JGR, 2007.
        !!       (3) BIO-LIM, Martin Vancoppenolle, 2008
        !!
        !!------------------------------------------------------------------
        !! * Arguments

      INCLUDE 'type.com'
      INCLUDE 'para.com'
      INCLUDE 'const.com'
      INCLUDE 'ice.com'
      INCLUDE 'thermo.com'
      INCLUDE 'bio.com'

      ! Energy conservation
      LOGICAL con_i

      zeps0 = 1.0e-16
      zeps1 = 1.0e-20
      zeps2 = 1.0e-04

      con_i = .true. ! conservation check in the ice or not
      jl = 1         ! category number (temporary)
      ji = 1         ! ji index

      WRITE(numout,*) ' * ice_th : '
      WRITE(numout,*) ' ~~~~~~~~~~ '
      WRITE(numout,*) 
      WRITE(numout,*) ' nlay_i : ', nlay_i
      WRITE(numout,*) ' nlay_s : ', nlay_s
      WRITE(numout,*) 
!
!-------------------------------------------------------------------------------
!  1) Initialize sea ice
!-------------------------------------------------------------------------------
!

      IF ( numit .EQ. nstart ) THEN
         CALL ice_phy_ini( 1 , 1 , nlay_s, nlay_i)

         IF ( c_bio_model .EQ. 'KRILL' )
     &   CALL ice_bio_ini( 1 , 1 , nlay_i )                     ! Bio initialization - KRILL

      ENDIF
!
!-------------------------------------------------------------------------------
!  2) Consistency & Energetic checks
!-------------------------------------------------------------------------------
!

      DO layer = 1, nlay_s
         t_s_b(ji,layer)  =  MIN( tpw , t_s_b(ji,layer) )
      END DO

      DO layer = 1, nlay_i
         tmelts           =  - tmut*s_i_b(ji,layer) + tpw
         t_i_b(ji,layer)  =  MIN(tmelts,t_i_b(ji,layer))
      END DO

      CALL ice_th_enmelt(1,1, nlay_s, nlay_i) ! ice enthalpy

      ! Initialize total heat content
      IF ( con_i ) CALL ice_th_glohec( qt_i_in , qt_s_in ,              
     &                                 q_i_layer_in , 1 , 1 , jl, 
     &                                 nlay_s , nlay_i )
!
!-------------------------------------------------------------------------------
!  3) Model routines
!-------------------------------------------------------------------------------
!
      CALL forcing                                              ! Compute forcing

      CALL ice_rad( nlay_s , nlay_i , 1 , 1 )                   ! Radiative transfer

      CALL ice_th_diff( nlay_s , nlay_i , 1 , 1 )               ! Heat diffusion

      IF ( con_i ) THEN                                         ! Conservation test
         CALL ice_th_glohec( qt_i_fin , qt_s_fin , q_i_layer_fin ,
     &                       1 , 1 , jl , nlay_s , nlay_i )
         CALL ice_th_con_dif( 1 , 1 , nlay_s , nlay_i , jl )
      ENDIF

      IF ( c_gravdr .EQ. 'CW' )
     &   CALL ice_sal_diff_CW(nlay_i,1,1)                       ! Salt transport (Cox and Weeks)

      IF ( c_gravdr .EQ. 'RA' )
     &   CALL ice_sal_diff(nlay_i,1,1)                          ! Salt transport (Rayleigh-number based)

      IF ( c_gravdr .EQ. 'AD' ) 
     &   CALL ice_sal_adv(nlay_i,1,1)                           ! Salt transport (Advection-based)

      IF ( ( c_bio_model .EQ. 'KRILL' ) .AND. ln_trdiff )
     &   CALL ice_bio_diff( 1 , 1 , nlay_i )                    ! Bio transport

      IF ( c_bio_model .EQ. 'KRILL' )
     &   CALL ice_bio_sms(nlay_i,1,1)                           ! KRILL source minus sinks
     
      IF ( ln_ikaite ) 
     &   CALL ice_ikaite(nlay_i)                                ! Ikaite precipitation and dissolution

      CALL ice_gas(nlay_i,1,1)                                  ! Gases

      CALL ice_th_dh(nlay_s,nlay_i,1,1)                         ! Growth and melt

      CALL ice_phy_remap(nlay_s,nlay_i,1,1)                     ! Remap heat and salt

      IF ( con_i ) THEN                                         ! Heat and salt conservation test
         CALL ice_th_glohec( qt_i_fin , qt_s_fin , q_i_layer_fin ,
     &                       1 , 1 , jl, nlay_s, nlay_i )
         CALL ice_th_con_dh(1,1,nlay_s,nlay_i,jl)
      ENDIF

      IF ( ln_trremp )      
     &   CALL ice_bio_remap(nlay_s, nlay_i, 1, 1)               ! Remap tracers

      CALL ice_bio_column(kideb,kiut,ntra_bio,ct_i_bio,cbu_i_bio,
     &                    deltaz_i_bio, .FALSE.)
      
      !---------------
      ! Chlorophyll a
      !---------------
      IF ( c_bio_model .EQ. 'KRILL' ) THEN
         DO layer = 1, nlay_bio
            chla_i_bio(layer) = cbu_i_bio(4,layer) * chlC_bio(layer) *
     &                          c_molar
         END DO
      ENDIF ! c_bio_model

      IF ( c_bio_model .NE. 'NOBIO' ) 
     &   WRITE(numout,*) ' chla_i_bio : ', ( chla_i_bio(layer), 
     &                   layer = 1, nlay_bio )

!-------------------------------------------------------------------------------
!  4) Outputs 
!-------------------------------------------------------------------------------
!
      !---------------
      ! Netcdf Output
      !---------------
      CALL ice_output(nlay_i,nlay_s)

      RETURN
!
!------------------------------------------------------------------------------
!- end of ice_th
      END