SUBROUTINE ice_bio_grid(kideb,kiut,nlay_i,ln_write)

! This routine creates biogeochemical vertical grid
! (c) Martin Vancoppenolle, May 2007
 
      INCLUDE 'type.com'
      INCLUDE 'para.com'
      INCLUDE 'const.com'
      INCLUDE 'ice.com'
      INCLUDE 'thermo.com'
      INCLUDE 'bio.com'

      INTEGER :: 
     &  ji          , ! : index for space
     &  jk          , ! : index for ice layers
     &  jn            ! : index for tracers

      LOGICAL ::
     &   ln_write

!=============================================================================!

      IF ( ln_write ) THEN
         WRITE(numout,*)
         WRITE(numout,*) ' *** ice_bio_grid : '
         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~ '
         WRITE(numout,*)
         WRITE(numout,*) ' c_grid : ', c_grid
      ENDIF

      DO ji = kideb, kiut

      !-----------------------------------------------------------------------!
      ! Thickness of the biologically-active layer: h_bio (SL & BA)
      !-----------------------------------------------------------------------!

      ! Skeletal layer case
      IF ( c_grid .EQ. 'SL' ) THEN
         h_bio = MIN( h_skel , ht_i_b(ji) )    !!! skeletal layer thickness is prescribed
      ENDIF

      ! Biologically-active layer case
      IF ( c_grid .EQ. 'BA' ) THEN

         ! Brine volume fraction
         DO layer = 1, nlay_i
            e_i_b(layer) = - tmut * s_i_b(ji,layer) / ( t_i_b(ji,layer) 
     &                      - tpw )
         END DO

         ! Find the layer above which BAL is located
         i_bal = 1
         DO layer = 2, nlay_i
            IF ( ( e_i_b(layer-1) .LT. e_thr_bal ) .AND. 
     &           ( e_i_b(layer)   .GT. e_thr_bal ) ) i_bal = layer
         END DO ! layer

         ! Upper-boundary of the BAL (linear interpolation)
         IF ( i_bal .GT. 1 ) THEN
            zdh = ht_i_b(ji) / REAL( nlay_i )
            zde = ( e_i_b( i_bal ) - e_i_b( i_bal-1 ) )
            zm  = zde / zdh
            zdz = ( e_thr_bal - e_i_b( i_bal - 1 ) ) / zm
            zbal = z_i_phy(i_bal-1) + zdz
         ELSE 
            zbal = 0.
         ENDIF

         ! BAL thickness
         h_bio = MIN( MAX( ht_i_b(ji) - zbal , h_skel ), 
     &                ht_i_b(ji) - h_skel ) !!! BAL layer is computed
         ! it cannot be higher than ht_i_b - h_skel
         ! it cannot be smaller than h_skel

      ENDIF ! c_grid .EQ. 'BA'

      !-----------------------------------------------------------------------!
      ! Layer thicknesses
      !-----------------------------------------------------------------------!

      ! Uniform multi-layer case
      IF ( c_grid .EQ. 'ML' ) THEN

         h_bio = ht_i_b(ji)                    !!! biologically-active layer = ice thickness

         zh = ht_i_b(ji) / REAL( nlay_bio )
         DO layer = 1, nlay_bio
            deltaz_i_bio(layer) = zh
         END DO

      ENDIF

      ! Deformed case (uppermost and lowermost have fixed thickness to h_skel)
      IF ( c_grid .EQ. 'DL' ) THEN
         h_bio = ht_i_b(ji)                    !!! biologically-active layer = ice thickness

         deltaz_i_bio(1)        = MIN( ht_i_b(ji) / REAL( nlay_bio ) , 
     &                            h_skel  )
         deltaz_i_bio(nlay_bio) = MIN( ht_i_b(ji) / REAL( nlay_bio ) ,
     &                            h_skel  )
         zh = ( ht_i_b(ji) - deltaz_i_bio(1) - deltaz_i_bio(nlay_bio) )
     &        / REAL( nlay_bio - 2 )
         DO layer = 2, nlay_bio - 1
            deltaz_i_bio(layer) = zh
         END DO
      ENDIF

      IF ( ( c_grid .EQ. 'SL' ) .OR. ( c_grid .EQ. 'BA' ) ) THEN
         deltaz_i_bio(nlay_bio) = h_bio
         DO layer = 1, nlay_bio - 1
            deltaz_i_bio(layer) = ( ht_i_b(ji) - h_bio ) / 
     &                            ( nlay_bio - 1 )
         END DO
         WRITE(numout,*) " h_bio : ", h_bio
         WRITE(numout,*) " ht_i_b(ji) : ", ht_i_b(ji)
         WRITE(numout,*) " deltaz_i_bio : ", ( deltaz_i_bio(layer), 
     &                   layer = 1, nlay_bio )
      ENDIF

      !-----------------------------------------------------------------------!
      ! Mid-point cotes
      !-----------------------------------------------------------------------!
      z_i_bio(1)  = deltaz_i_bio(1) / 2.0
      DO layer = 2, nlay_bio
         z_i_bio(layer)  = z_i_bio(layer-1) + ( deltaz_i_bio(layer-1) +
     &                     deltaz_i_bio(layer) ) / 2.0
      END DO

      !-----------------------------------------------------------------------!
      ! Layer boundaries
      !-----------------------------------------------------------------------!
      zb_i_bio(0) = 0.0
      DO layer = 1, nlay_bio
         zb_i_bio(layer) = zb_i_bio(layer-1) + deltaz_i_bio(layer)
      END DO ! layer

      !-----------------------------------------------------------------------!
      END DO ! ji

      IF ( ln_write ) THEN
         WRITE(numout,*) ' deltaz_i_bio : ', ( deltaz_i_bio(layer),  
     &                   layer = 1, nlay_bio)
         WRITE(numout,*) ' z_i_bio  : ', ( z_i_bio(layer),  
     &                   layer = 1, nlay_bio)
         WRITE(numout,*) ' zb_i_bio : ', ( zb_i_bio(layer),  
     &                   layer = 1, nlay_bio)

         WRITE(numout,*)
      ENDIF

!=============================================================================!
!-- End of ice_bio_ini --
 
      END