SUBROUTINE ice_ikaite(nlay_i)

!------------------------------------------------------------------------------!
!                               *** ice_ikaite ***
!
!     This routine computes precipitation / dissolution of ikaite 
!     and impacts on DIC, Alk and Ca2+
!
!     Ref: Moreau et al (JGR2015)
!     
!     Original code: S. Moreau, M. Vancoppenolle, 2012-2015
!     Rewriting: M. Vancopppenolle, Nov 2015
!
!     Solubility product formulation is valid until -9C
!     The fit is not monotonic below -9C, so we may experience problems 
!     with ikaite below -9C
!
!------------------------------------------------------------------------------!

      INCLUDE 'type.com'
      INCLUDE 'para.com'
      INCLUDE 'const.com'
      INCLUDE 'ice.com'
      INCLUDE 'thermo.com'
      INCLUDE 'bio.com'
      
       REAL(8), DIMENSION(nlay_bio) :: 
     &    Kspi

       LOGICAL ::
     &   ln_write_bio = .TRUE.

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

!
!---------------------------------------------------------------------
! 1) Starting the routine
!---------------------------------------------------------------------
!

      IF ( ln_write_bio ) THEN

         WRITE(numout,*)
         WRITE(numout,*) ' *** ice_ikaite : '
         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~'
         WRITE(numout,*)
         WRITE(numout,*) ' *** Initial values before CaCO3 precip/diss '
         WRITE(numout,*) '        at time step : ', numit

      ENDIF
      
      CALL ice_brine     ! get brine salinity and density
      
      CALL ice_carb_chem ! get CO32m 

      DO jn = 1, ntra_bio

         ! Brine / Total gas concentrations
         DO layer = 1, nlay_bio
             c_i_bio(jn,layer)   = cbu_i_bio(jn,layer) / e_i_bio(layer)
         END DO

         IF ( ln_write_bio ) THEN

            WRITE(numout,*)
            WRITE(numout,*) ' Tracer :   ', biotr_i_nam(jn)
            WRITE(numout,*) ' c_i_bio: ', ( c_i_bio(jn,layer),  ! concentration of tracers in bulk ice
     &                                    layer = 1, nlay_bio )
            WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer),  ! concentration of tracers in bulk ice
     &                                      layer = 1, nlay_bio )
            WRITE(numout,*) 
         ENDIF
      END DO
      
!
!------------------------------------------------------------------------------!
! 2) Solubility product
!------------------------------------------------------------------------------!
!
      ! Ref: Papadimitriou et al. (GCA 2013)
      ! fit valid until -9 degC
      ! non monotonic (at all) below this value
      ! we may have to regularize this
      
      DO layer = 1, nlay_bio

          zT  = t_i_bio(layer)

          IF ( zt .LE. 273.15-9.) zT = 273.15-9. ! asked stathys papadim what to do here

          zpKsp = - 15489.09608 + 623443.70216 / zT +
     &             2355.14596 * log(zT)

          Kspi(layer) = 10.**(-zpKsp)     ! (mol/kg)^2 of brine
          
          ! Conversion into (mmol/m3)^2
          Kspi(layer) = Kspi(layer) * 1.0e6 * 
     &                  rhobr_bio(layer) * rhobr_bio(layer)

      END DO
        
      WRITE(numout,*) ' Kspi : ', Kspi(1:nlay_bio)
      
!
!------------------------------------------------------------------------------!
! 3) Calcification / Dissolution
!------------------------------------------------------------------------------!
!

      DO layer = 1, nlay_bio

         zcal = c_i_bio(jn_cal,layer)
         ! zcal = 311.765 * sbr_bio(layer) --> approach not valid because Ca2+ changes a lot in brine in winter

         ! Brine concentration of CO32-
         zco32m_br  = co32m(layer) / e_i_bio(layer)
         zco32m_sat = Kspi(layer)  / zcal
         
         ! Saturation state (omega) in brine
         zomega = zco32m_br * zcal / Kspi(layer)
         
         ! Precipitation / dissolution
         zd_ika = ddtb / caco3_time * ( zco32m_br - zco32m_sat )
         zd_ika = MAX( MIN(zd_ika, zco32m_br), -c_i_bio(jn_ika,layer) ) ! calcification cannot exceed CO32- stock
                                                                        ! dissolution cannot exceed ikaite stock 
         zd_ika = zd_ika * e_i_bio(layer)                               ! bulk calcification 

         ! Update ikaite, dic, alk and calcium bulk concentrations
         cbu_i_bio(jn_ika,layer) = cbu_i_bio(jn_ika,layer) + zd_ika
         cbu_i_bio(jn_dic,layer) = cbu_i_bio(jn_dic,layer) - zd_ika
         cbu_i_bio(jn_alk,layer) = cbu_i_bio(jn_alk,layer) - 2.*zd_ika
         cbu_i_bio(jn_cal,layer) = cbu_i_bio(jn_cal,layer) - zd_ika

         ! Retrieve brine concentrations
         c_i_bio(jn_ika,layer) = cbu_i_bio(jn_ika,layer)/e_i_bio(layer)
         c_i_bio(jn_dic,layer) = cbu_i_bio(jn_dic,layer)/e_i_bio(layer)
         c_i_bio(jn_alk,layer) = cbu_i_bio(jn_alk,layer)/e_i_bio(layer)
         c_i_bio(jn_cal,layer) = cbu_i_bio(jn_cal,layer)/e_i_bio(layer)

         ! diagnostic rate of formation and saturation state
         ika_rate(layer)  = zd_ika / ddtb
         ika_omega(layer) = zomega

         WRITE(numout,*) ' layer, zd_ika : ', layer, zd_ika

      END DO
      
!
!------------------------------------------------------------------------------!
! 5) End of the routine
!------------------------------------------------------------------------------!
!

      WRITE(numout,*) 
      WRITE(numout,*) ' ** After CaCO3 precipitation : '
      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
      WRITE(numout,*)

       DO jn = 1, ntra_bio

         IF ( ( jn .EQ. jn_dic ) .OR.
     &        ( jn .EQ. jn_alk ) .OR.
     &        ( jn .EQ. jn_cal ) .OR.
     &        ( jn .EQ. jn_ika ) ) THEN
            WRITE(numout,*) ' Tracer :   ', biotr_i_nam(jn)
            WRITE(numout,*) ' c_i_bio  : ', ( c_i_bio(jn,layer),    ! concentration of tracers in brines (mmol m-3)
     &                                        layer = 1, nlay_bio )
            WRITE(numout,*) ' cbu_i_bio: ', ( cbu_i_bio(jn,layer),  ! concentration of tracers in bulk ice
     &                                        layer = 1, nlay_bio )
         ENDIF

      END DO

      WRITE(numout,*)
      WRITE(numout,*) ' End of ice_ikaite '
      WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
     

!==============================================================================|
! End of ice_ikaite

      WRITE(numout,*)
      END