SUBROUTINE ice_bio_remap(nlay_s,nlay_i,kideb,kiut)
!-----------------------------------------------------------------------------!
      INCLUDE 'type.com'
      INCLUDE 'para.com'
      INCLUDE 'const.com'
      INCLUDE 'ice.com'
      INCLUDE 'thermo.com'
      INCLUDE 'bio.com'

      REAL(8), DIMENSION( 0:maxnlay+2 ) ::
     &  z0

      REAL(8), DIMENSION( maxnlay + 2 ) ::
     &  zq0         , ! : scalar content on the physical grid (input)
     &  zthick0       ! : thickness of biological layers

      REAL(8), DIMENSION( maxnlay ) ::
     &  zq1         , ! : scalar content on the biological grid (output)
     &  zdeltaz

      REAL(8), DIMENSION( nlay_bio , maxnlay ) ::
     &  zweight       ! : relayering matrix

      INTEGER :: zcase
      zeps   = 1.0e-20

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

      WRITE(numout,*) 
      WRITE(numout,*) ' ice_bio_remap : '
      WRITE(numout,*) ' ~~~~~~~~~~~~~~~ '
      WRITE(numout,*)
      WRITE(numout,*) '    *** Before remapping                   *** '
      DO jn = 1, ntra_bio
         IF ( flag_active(jn) ) THEN
           WRITE(numout,*) ' biotr_i_nam : ', biotr_i_nam(jn)
           WRITE(numout,*) ' cbu_i_bio : ', ( cbu_i_bio(jn, jk), jk = 1,
     &                                      nlay_bio )
         ENDIF
      END DO

      !--------------------
      ! Conservation check
      !--------------------
      CALL ice_bio_column(kideb,kiut,ntra_bio,mt_i_bio_init,cbu_i_bio,
     &                    deltaz_i_bio, .FALSE.)

      DO ji = kideb, kiut
      WRITE(numout,*)
      WRITE(numout,*) ' Inputs     : '
      WRITE(numout,*) ' ~~~~~~       '
      WRITE(numout,*) ' dh_i_surf  : ', dh_i_surf(ji)
      WRITE(numout,*) ' dh_i_bott  : ', dh_i_bott(ji)
      WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji)
!
!------------------------------------------------------------------------------|
!  1) Switches                                                                 |
!------------------------------------------------------------------------------|
!
      WRITE(numout,*)
      WRITE(numout,*) ' Switches : '
      WRITE(numout,*) ' ~~~~~~~~ '

      !----------------------------------
      ! Biological Grid Type
      !----------------------------------
      IF ( ( c_grid .NE. 'SL' ) .AND. ( c_grid .NE. 'BA' ) ) THEN
         zcase = 1 ! ML / DL cases
      ELSE
         zcase = 2 ! SL / BA cases
      ENDIF

      !----------------------------------
      ! Surface melt indicator : icsuind
      !----------------------------------
      ! ice surface behaviour : computation of icsuind-icsuswi
      ! icsuind : index which values equals 
      !           0 if there is nothing
      !           1 if first layer has started to melt
      !           2 if first and second layer have started ...
      !                             etc ...
      icsuind  =  0  
      deltah   =  0.0
      DO layer = 1, nlay_bio
         IF ( - dh_i_surf(ji) .GT. deltah ) THEN
            icsuind  =  layer 
         ENDIF
         deltah = deltah + deltaz_i_bio(layer)
      END DO

      !-------------------------------
      ! Surface melt switch : icsuswi
      !-------------------------------
      ! icsuswi : switch which value equals 1 if ice melts at the surface
      !                                     0 if not
      icsuswi    =  MAX( 0 , INT( - dh_i_surf(ji) 
     &              / MAX ( zeps , ABS( dh_i_surf(ji) ) ) ) )

      !--------------------------------
      ! Ice bottom indicator : icboind
      !--------------------------------
      ! icboind : index which values equals 0 if accretion is on the way
      !                                     1 if last layer has started to melt
      !                                     2 if penultiem layer has started ...
      !                                     etc ...
      icboind  =  0  
      deltah   =  0.0
      DO layer = nlay_bio, 1, -1
         IF ( - dh_i_bott(ji) .GT. deltah ) THEN
            icboind  =  nlay_bio + 1 - layer 
         ENDIF
         deltah = deltah + deltaz_i_bio(layer)
      END DO

      !-----------------------------
      ! Ice bottom switch : icboswi
      !-----------------------------
      ! icboswi : switch which value equals 1 if accretion of ice is on the way
      !                                     0 if ablation is on the way
      icboswi    =  MAX( 0 , INT( dh_i_bott(ji)
     &              / MAX( zeps , ABS( dh_i_bott(ji) ) ) ) )

      !------------------------------
      ! Snow ice indicator : snicind
      !------------------------------
      ! snow-ice formation : calcul de snicind-snicswi
      ! snicind : index which values equals 0 if no snow-ice forms
      !                1 if last layer of snow has started to melt
      !                2 if penultiem layer ...
      snicind  =  0
      deltah   =  0.0
      DO layer = nlay_s, 1, -1
         IF ( dh_snowice(ji) .GT. deltah) THEN
            snicind  =  nlay_s+1-layer
         ENDIF
         deltah = deltah + h_s(layer)
      END DO

      !---------------------------
      ! Snow ice switch : snicswi
      !---------------------------
      ! snicswi : switch which value equals 1 if snow-ice forms
      !                                     0 if not
      snicswi    =  MAX( 0 , INT( dh_snowice(ji) /
     &              MAX( zeps , ABS( dh_snowice(ji) ) ) ) )

      WRITE(numout,*) ' icsuind    : ', icsuind
      WRITE(numout,*) ' icsuswi    : ', icsuswi
      WRITE(numout,*) ' icboind    : ', icboind
      WRITE(numout,*) ' icboswi    : ', icboswi
      WRITE(numout,*) ' snicind    : ', snicind
      WRITE(numout,*) ' snicswi    : ', snicswi
      WRITE(numout,*) 
!
!------------------------------------------------------------------------------|
! 2) Tracer sources
!------------------------------------------------------------------------------|
!
      WRITE(numout,*) ' Tracer sources : '
      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
      
      DO jn = 1, ntra_bio

         IF ( flag_active (jn) ) THEN

         c_s_bio(jn)   = 0.0 ! now no concentration in the snow
         ch_bo_bio(jn) = 0.0
         ch_si_bio(jn) = 0.0

         IF ( ln_trbo ) 
     &      ch_bo_bio(jn) = e_skel * c_w_bio(jn) * MAX( dh_i_bott(ji),
     &                      0.0 )

         IF ( ln_trsi )
     &      ch_si_bio(jn) = ( ( rhog - rhon ) / rhog * c_w_bio(jn) 
     &                   + ( rhon / rhog ) * c_s_bio(jn) ) * 
     &                   MAX( dh_snowice(ji) , 0.0 )

         fcbp(jn) = ch_bo_bio(jn) / ddtb
         fcsi(jn) = ch_si_bio(jn) / ddtb

         WRITE(numout,*) ' Tracer    : ', biotr_i_nam(jn)
         WRITE(numout,*) ' ch_bo_bio : ', ch_bo_bio(jn)
         WRITE(numout,*) ' ch_si_bio : ', ch_si_bio(jn)
!        WRITE(numout,*) ' c_s_bio   : ', c_s_bio (jn)
!        WRITE(numout,*) ' c_w_bio: ', c_w_bio(jn)
!        WRITE(numout,*) ' c_ws_i    : ', 
!    &                   ( rhog - rhon ) / rhog * c_w_bio(jn)

         ENDIF

      END DO !jn
!
!------------------------------------------------------------------------------|
! 3) Compute new concentrations by remapping tracers
!------------------------------------------------------------------------------|
!
      SELECT CASE (zcase) 
                
                 !========================!
      CASE (1)   !=== MULTI-LAYER CASE ===!
                 !========================!

      !
      !------------------------------------------------------------------------|
      ! 3.1) Old grid: zthick0, z0, zdeltaz
      !------------------------------------------------------------------------|
      !
      WRITE(numout,*) ' Old grid : '
      WRITE(numout,*) ' ~~~~~~~~~~ '

      !--- Diffusive remapping, compute z0 and zthick0 ---!
      !
      ! indexes of the vectors
      ntop0    =  1
      nbot0    =  nlay_bio + 1 - icboind + (1 - icsuind) * icsuswi +
     &            snicswi
      nbot0    =  MAX(1, MIN( nlay_bio + 1 - icboind + 
     &            (1-icsuind)*icsuswi
     &          + snicswi, nlay_bio + 2 ) )

      WRITE(numout,*) ' ntop0     : ', ntop0
      WRITE(numout,*) ' nbot0     : ', nbot0

      ! Cotes of the top of the layers
      z0(0)   =  0.0
      DO layer = 1, nbot0-1
        limsum    =  ( (icsuswi * (icsuind+layer-1) + 
     &                 (1-icsuswi)*layer) ) *
     &                 (1-snicswi) + (layer-1)*snicswi
         z0(layer) =  icsuswi*dh_i_surf(ji) + snicswi*dh_snowice(ji)
         DO layer_a = 1, limsum
            z0(layer) = z0(layer) + deltaz_i_bio(layer_a)
         END DO
      END DO

      z0(nbot0)  =  icsuswi*dh_i_surf(ji) + snicswi*dh_snowice(ji) + 
     &              dh_i_bott(ji)

      DO layer = 1, nlay_bio
         z0(nbot0) = z0(nbot0) + deltaz_i_bio(layer)
      END DO
       
      z0(1)      =  snicswi * dh_snowice(ji) + (1 - snicswi) * z0(1)

      DO layer = ntop0, nbot0
         zthick0(layer)  =  z0(layer) - z0(layer-1)
      END DO

      WRITE(numout,*) ' Old grid, ntop0 : ', ntop0, ' nbot0 : ', nbot0
      WRITE(numout,*) ' z0      : ', ( z0(layer), layer = 0, nbot0 )
      WRITE(numout,*) ' zthick0 : ', ( zthick0(layer) ,
     &                               layer = 1, nbot0 )
      WRITE(numout,*)

      !--- Squeezing remapping, compute zdeltaz ---!'
      !
      zdeltaz(1:nlay_bio) = deltaz_i_bio(:)

      !
      !-------------------------------------------------------------------------|
      ! 3.2) New grid   
      !-------------------------------------------------------------------------|
      !
      WRITE(numout,*) ' New grid : '
      WRITE(numout,*) ' ~~~~~~~~~~ '
      ! indexes of the vectors
      ntop1    =  1 
      nbot1    =  nlay_bio
      ! cotes and thicknesses
      CALL ice_bio_grid(kideb,kiut,nlay_i,.TRUE.) ! compute the biological grid

      WRITE(numout,*) ' New grid, ntop1 : ', ntop1, ' nbot1 : ', nbot1
      WRITE(numout,*) ' zb_i_bio: ', ( zb_i_bio(layer), 
     &                                 layer = 0, nbot1 )
      !
      !------------------------------------------------------------------------|
      ! 3.3) Weights (diffusive remapping only)
      !------------------------------------------------------------------------|
      !
      WRITE(numout,*)
      WRITE(numout,*) ' Weights  : '
      WRITE(numout,*) ' ~~~~~~~~~~ '

      DO layer1 = 1, nlay_bio
         DO layer0 = 1, nbot0
            zweight(layer1,layer0) = MAX ( 0.0 , ( MIN ( z0(layer0) ,
     &      zb_i_bio(layer1) ) 
     &    - MAX ( z0 (layer0-1) , zb_i_bio(layer1-1) ) ) / 
     &      zthick0(layer0) )
!           WRITE(numout,*) ' layer0, layer1 : ', layer0, layer1
!           WRITE(numout,*) ' zweight        : ', zweight(layer1,layer0)
         END DO
      END DO
      WRITE(numout,*)

      !
      !------------------------------------------------------------------------|
      ! 3.4) Tracer remapping
      !------------------------------------------------------------------------|
      !
      WRITE(numout,*) ' Tracer remapping : '
      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~ '

      DO jn = 1, ntra_bio

      IF ( flag_active (jn) ) THEN

         WRITE(numout,*) ' ------------------------------------------'
         WRITE(numout,*) ' Tracer                  : ', biotr_i_nam(jn)
         WRITE(numout,*) ' Remapping type  nn_remp : ', nn_remp(jn)
         WRITE(numout,*) ' ------------------------------------------'
         WRITE(numout,*)

         !--------------------------------------------------------------------!
         ! 3.4.1) Compute 'new' tracer contents
         !--------------------------------------------------------------------!

                                          !=====================!
         IF ( nn_remp(jn) .EQ. 1 ) THEN   ! Diffusive remapping !
                                          !=====================!
         !--------------------!
         ! Old tracer content !
         !--------------------!

         WRITE(numout,*) ' Tracer content ...'
         WRITE(numout,*)

         !--- Inner layers ---!
         mt_i_bio_init(jn) = 0.0

         DO layer = ntop0, nbot0
            limsum =  MAX(1,MIN(snicswi*(layer-1) + icsuswi*(layer-1
     &           +  icsuind)
     &           + (1-icsuswi)*(1-snicswi)*layer, nlay_bio))
            zq0(layer) = zthick0(layer) * cbu_i_bio(jn,limsum)

!           WRITE(numout,*) ' limsum : ', limsum
!           IF ( limsum .GE. 1 ) 
!    &         WRITE(numout,*) ' cbu_i_bio: ', cbu_i_bio(jn,limsum)

!           WRITE(numout,*) ' zq0    : ', zq0(layer)
!           WRITE(numout,*) ' zthick0: ', zthick0(layer)

         END DO

         !--- Bottom layer ---!
         ! bottom layer if ice forms at the bottom
         zq0(nbot0) = ( 1 - icboswi ) * zq0(nbot0) + icboswi *
     &                ch_bo_bio(jn)

         !--- Upper layer ---!
         ! first layer if snow ice forms
         zq0(1)     =  snicswi * ch_si_bio(jn) 
     &              + ( 1 - snicswi ) * zq0(1)

!        WRITE(numout,*) ' snicswi   : ', snicswi
!        WRITE(numout,*) ' c_si_bio  : ', ch_si_bio(jn) / zthick0(1)

         DO layer = ntop0, nbot0
            mt_i_bio_init(jn) = mt_i_bio_init(jn) +
     &                          zq0(layer)
         END DO

         WRITE(numout,*) ' zq0       : ', ( zq0(layer), 
     &                                    layer = 1, nbot0 )
         !----------------!
         ! Remapping      !
         !----------------!

         WRITE(numout,*) ' Remapping ... ' 
         WRITE(numout,*)

         DO layer1 = 1, nlay_bio
            zq1(layer1) = 0.0
            DO layer0 = 1, nbot0
               zq1(layer1) = zq1(layer1) + zweight(layer1,layer0) *
     &                       zq0(layer0)
            END DO
         END DO

         ENDIF ! nn_remp

                                          !=====================!
         IF ( nn_remp(jn) .EQ. 2 ) THEN   ! Squeezing remapping
                                          !=====================!
 
         WRITE(numout,*) ' Remapping ... ' 
         WRITE(numout,*)


         ! Inner layers
         DO layer = 1, nlay_bio 
            zq0(layer) = zdeltaz(layer) * cbu_i_bio(jn,layer)
         END DO

         ! Surface layer
         zq0(1)        = zq0(1) + ch_si_bio(jn)

         ! Bottom layer
         zq0(nlay_bio) = zq0(nlay_bio) + ch_bo_bio(jn)

         ! Tracer content for conservation check
         mt_i_bio_init(jn) = 0.0
         DO layer = 1, nlay_bio
            mt_i_bio_init(jn) = mt_i_bio_init(jn) +
     &                          zq0(layer)
         END DO

         zq1(1:nlay_bio) = zq0(1:nlay_bio)

         ENDIF ! nn_remp

         WRITE(numout,*)
         WRITE(numout,*) ' zq1       : ', ( zq1(layer), 
     &                                    layer = 1, nlay_bio )
         WRITE(numout,*)

         !--------------------------------------------------------------------!
         ! 3.4.2) Retrieve bulk concentration from content
         !--------------------------------------------------------------------!
         !
         DO layer = 1, nlay_bio
            cbu_i_bio(jn,layer) = zq1(layer) / deltaz_i_bio(layer)
         END DO ! layer
      
      ENDIF ! flag_active

      END DO ! jn

                 !========================!
      CASE (2)   !=== SL/BAL CASE         !
                 !========================!

         zdz0 = deltaz_i_bio(nlay_bio)

         CALL ice_bio_grid(kideb,kiut,nlay_i,.TRUE.) ! compute the new biological grid

         DO jn = 1, ntra_bio
         
            IF ( flag_active (jn) ) THEN
      
            ! Old tracer content
            zq0(1:nlay_bio-1) = 0.
            zq0(nlay_bio) = snicswi * ch_si_bio(jn) 
     &                    + icboswi * ch_bo_bio(jn)
     &                    + cbu_i_bio(jn,nlay_bio) * zdz0

            ! New tracer content (no loss from melt)
            zq1(1:nlay_bio-1) = 0.
            zq1(nlay_bio)     = zq0(nlay_bio)

            ! New bulk concentrations
            cbu_i_bio(jn,1:nlay_bio-1) = 0.0
            cbu_i_bio(jn,nlay_bio)     = zq1(nlay_bio) / 
     &                                   deltaz_i_bio(nlay_bio)

            ENDIF ! flag_active

         END DO ! jn

      END SELECT ! zcase

!
!-----------------------------------------------------------------------
! 4) Retrieve brine concentration
!-----------------------------------------------------------------------
!
      !---------------
      ! Interpolation 
      !---------------
      CALL ice_bio_interp_phy2bio(kideb,kiut,nlay_i,.FALSE.) 
                                          ! interpolation of physical variables
                                          ! on the biological grid
      !---------------------
      ! Brine concentration
      !---------------------
      DO jn = 1, ntra_bio

         IF ( flag_active (jn) .AND. flag_diff(jn) ) THEN
            DO layer = 1, nlay_bio
               c_i_bio(jn,layer) = cbu_i_bio(jn,layer) / e_i_bio(layer)
            END DO ! layer
         ENDIF ! flags

      END DO ! jn

      END DO ! ji

!
!-----------------------------------------------------------------------
! 5) Conservation check
!-----------------------------------------------------------------------
!
      !-------------------
      ! Conservation test
      !-------------------
      CALL ice_bio_column(kideb,kiut,ntra_bio,mt_i_bio_final,cbu_i_bio,
     &                    deltaz_i_bio, .FALSE.)
      DO jn = 1, ntra_bio
         IF ( flag_active(jn) ) THEN
            WRITE(numout,*) ' mt_i_bio_final : ', mt_i_bio_final(jn)
         ENDIF
      END DO

      DO jn = 1, ntra_bio

         IF ( flag_active(jn) ) THEN

         ! Flux due to bottom formation
         f_bogr(jn) =  ch_bo_bio(jn) / ddtb
         f_sigr(jn) =  ch_si_bio(jn) / ddtb
         f_bogr(jn) =  0.0
         f_sigr(jn) =  0.0

         ! Flux due to snow ice formation

         ! Bottom flux ( positive upwards )
!        f_bo_tra(jn) =
!        f_bo_tra(jn) = zswitchw * ( - e_i_bio( nlay_bio )
!    &                * diff_bio * 2.0
!    &                / deltaz_i_bio(nlay_bio) * ( c_i_bio(jn,nlay_bio)
!    &                  - c_w_bio(jn) ) )
!    &                + zswitchs * ( - qsummer * c_i_bio(jn,nlay_bio) )
!    &                / ddtb
!    &                * e_i_bio(nlay_bio)
!        f_su_tra(jn) = zswitchw * 0.0
!    &                + zswitchs * ( qsummer * c_s_bio(jn) )
!    &                / ddtb
         !+++++
         WRITE(numout,*) ' f_bogr : ', f_bogr(jn)
         WRITE(numout,*) ' f_sigr : ', f_sigr(jn)
         !+++++

         ENDIF

      END DO ! jn

      zerror = 1.0e-9
      CALL ice_bio_conserv(kideb,kiut,ntra_bio,'ice_bio_remap  ',
     &                     biotr_i_nam, zerror,
     &                           mt_i_bio_init,mt_i_bio_final,
     &                           f_bogr, f_sigr, ddtb)
!     CALL ice_bio_column(kideb,kiut,ntra_bio,ct_i_bio,cbu_i_bio,
!    &                    deltaz_i_bio, .FALSE.)

!
!-----------------------------------------------------------------------
! 9) Final output
!-----------------------------------------------------------------------
!

      WRITE(numout,*)
      WRITE(numout,*) '    *** After remapping                    *** '
      WRITE(numout,*) '    model output '
      DO jn = 1, ntra_bio
         IF ( flag_active(jn) ) THEN
           WRITE(numout,*) ' biotr_i_nam : ', biotr_i_nam(jn)
           WRITE(numout,*) ' cbu_i_bio : ', ( cbu_i_bio(jn,jk),
     &                                      jk = 1, nlay_bio )
           IF ( flag_diff(jn) ) 
     &     WRITE(numout,*) ' c_i_bio   : ', ( c_i_bio(jn,jk),
     &                                      jk = 1, nlay_bio )
         ENDIF
      END DO

      WRITE(numout,*)
      WRITE(numout,*) ' End of ice_bio_remap '

!==============================================================================|
! Fin de la routine ice_bio_remap
      WRITE(numout,*)

      END