SUBROUTINE ice_gas(nlay_i,kideb,kiut)

!
!-----------------------------------------------------------------------------!
!                           *** ice_gas ***
!     
!     Nov 2015
!
!     Version : source_3.10
!
!     to do
!     - put the computation of solubilities in a subroutine
!        - check conservation
!        - brine characteristics to subroutinize    
!        - missing carbonate chemistry in the end
!-----------------------------------------------------------------------------!

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

       REAL(8) :: 
     &   zsat,
     &   zsat_diff

       LOGICAL ::
     &   ln_write_gas = .TRUE.

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

      IF ( ln_write_gas ) THEN

         WRITE(numout,*) 
         WRITE(numout,*) ' ** ice_gas : '
         WRITE(numout,*) ' ~~~~~~~~~~~~ '
   
      ENDIF
      
!
!------------------------------------------------------------------------------!
! 1) Initialize conservation test and gas concentrations
!------------------------------------------------------------------------------!
!

      !--- Gas concentrations
      DO jn = 1, ntra_bio

         ! Brine concentrations
         IF ( ( flag_active(jn)            ) .AND.
     &        ( biotr_i_typ(jn) .EQ. 'gas' ) ) THEN
            DO layer = 1, nlay_bio
              c_i_bio(jn,layer)   = cbu_i_bio(jn,layer) / e_i_bio(layer)
            END DO
         ENDIF

         ! Carbonate chemistry
         IF ( ln_carbon ) CALL ice_carb_chem

         IF ( ( flag_active(jn)            ) .AND.
     &        ( biotr_i_typ(jn) .EQ. 'gas' ) ) THEN
            DO layer = 1, nlay_bio
              c_gtot_i(jn,layer)  = cbu_i_bio(jn,layer) + 
     &                              cbub_i_bio(jn,layer)
            END DO

         ENDIF

      END DO

      !--- Conservation test
      CALL ice_gas_column( ntra_bio , c_gtot_i, deltaz_i_bio, .TRUE.,  
     &                     mt_i_gas_init )

      DO jn = 1, ntra_bio

         IF ( ( flag_active(jn)            ) .AND.
     &        ( biotr_i_typ(jn) .EQ. 'gas' ) .AND. 
     &        ln_write_gas                   ) THEN

           WRITE(numout,*)
           WRITE(numout,*) ' *** Initial values  *** '
           WRITE(numout,*) ' --- time step --- : ', numit

           WRITE(numout,*) ' --- Tracer    --- :   ', biotr_i_nam(jn)
           WRITE(numout,*) 'c_i_bio  : ',c_i_bio(jn,layer_00:nlay_bio)
           WRITE(numout,*) 'cbu_i_bio: ',cbu_i_bio(jn,layer_00:nlay_bio)
           WRITE(numout,*) 'cbub_i_bio: ',
     &                     cbub_i_bio(jn,layer_00:nlay_bio)
           WRITE(numout,*) 'c_gtot_i : ',c_gtot_i(jn,layer_00:nlay_bio)
           WRITE(numout,*)
           WRITE(numout,*) ' mt_i_gas_init : ', mt_i_gas_init(jn)

         ENDIF

      END DO
      
!
!------------------------------------------------------------------------------!
! 2) Bubble formation and dissolution
!------------------------------------------------------------------------------!
!

                          !---------------------------
      CALL ice_gas_solu   ! Saturation concentrations
                          !---------------------------
                  
      IF ( ln_bubform ) THEN
      
         IF ( ln_write_gas ) THEN
            WRITE(numout,*)
            WRITE(numout,*) ' *** Bubble formation and dissolution *** '
            WRITE(numout,*)
         ENDIF
         
          DO jn = 1, ntra_bio
          
             IF ( flag_active(jn) .AND. 
     &         ( biotr_i_typ(jn) .EQ. 'gas' ) ) THEN
      
                IF ( ln_write_gas ) THEN
                   WRITE(numout,*) ' Tracer number : ', jn
                   WRITE(numout,*) ' --- ', biotr_i_nam(jn)
                ENDIF
      
                DO layer = 1, nlay_bio
                
                   zsat      = csat_gas(jn,layer) * sursat_gas
                   zsat_diff = zsat - c_i_bio(jn,layer)    ! undersaturation
                   
                   ! bubble formation or dissolution rate
                   IF ( zsat_diff .LE. 0 ) THEN  ! sursaturation
                      zbub_rate = bub_form_rate
                   ELSE
                      zbub_rate = bub_diss_rate
                   ENDIF
                   
                   ! change in bulk bubble concentration
                   zd_cbri = zsat_diff * zbub_rate * ddtb
                   zd_cbri = MAX( -c_i_bio(jn,layer), zd_cbri )         ! cannot exhaust brines
                   zd_cbri = MIN( cbub_i_bio(jn,layer)/e_i_bio(layer),  ! cannot exhaust gas bubbles
     &                            zd_cbri )
                   c_i_bio(jn,layer) = c_i_bio(jn,layer) + zd_cbri
                
                   ! change in bulk bubble concentration
                   zd_cbub = - zd_cbri*e_i_bio(layer)
                   cbub_i_bio(jn,layer) = cbub_i_bio(jn,layer) + zd_cbub
 
                   ! Update bulk concentration
                   cbu_i_bio(jn,layer) =e_i_bio(layer)*c_i_bio(jn,layer)
 
                   ! Diagnostic % saturation ratio
                   pc_sat(jn,layer) = c_i_bio(jn,layer) / 
     &                                csat_gas(jn,layer) * 100.
 
                   ! Update DIC for CO2
                   IF ( jn .EQ. jn_co2 ) 
     &               cbu_i_bio(jn_dic,layer) = 
     &                  cbu_i_bio(jn_dic,layer) - zd_cbub
             
                END DO ! layer
          
             ENDIF ! flags, biotr_i_typ
             
          END DO ! ntra_bio
!
       ENDIF ! ln_bubform

!
!------------------------------------------------------------------------------!
! 3) Bubble rise
!------------------------------------------------------------------------------!
!
      IF ( ln_bubrise ) THEN
      
         IF ( ln_write_gas ) THEN
            WRITE(numout,*)
            WRITE(numout,*) ' *** Bubble rise *** '
         ENDIF
         
         f_bub(:) = 0.
          
         IF ( ln_write_gas ) WRITE(numout,*) ' *** immediate rise'
         
             DO jn = 1, ntra_bio
         
                IF ( flag_active(jn) .AND. 
     &             ( biotr_i_typ(jn) .EQ. 'gas' ) ) THEN
     
                   IF ( ln_write_gas ) THEN
                      WRITE(numout,*) ' Tracer number : ', jn
                      WRITE(numout,*) ' --- ', biotr_i_nam(jn)
                   ENDIF
                   
                   DO layer = nlay_bio, 1, -1
                   
                   zswitch = 0. ! by default, bubbles rise
                   
                   ! i_imper_layer: identifies the first impermeable layer above "layer"
                   ! if no such layer, then i_imper_layer = 0 and gas get out of the ice
                   ! via f_bub
                   
                   i_imper_layer = 0
                   
                   IF ( e_i_bio(layer) .LE. e_thr_bubrise ) THEN 
                      i_imper_layer = layer
                   ELSE
                      DO layer_2 = layer - 1, 1, -1 ! find a layer above
                       IF ( ( zswitch .LE. 0 ) .AND. 
     &                    ( e_i_bio(layer_2) .LT. e_thr_bubrise ) ) THEN
                            zswitch = 1.
                            i_imper_layer = MAX(i_imper_layer, layer_2)
                       ENDIF ! e_i_bio
                      END DO ! layer2
                   ENDIF

                   ! Transfer gas among layers
                   IF ( ( i_imper_layer .GT. 0 ) .AND. 
     &                  ( i_imper_layer .NE. layer ) ) THEN ! in-ice bubble transfer
                      cbub_i_bio(jn,i_imper_layer) = ! add gas in the impermeable layer
     &                cbub_i_bio(jn,i_imper_layer)+cbub_i_bio(jn,layer)
                      cbub_i_bio(jn,layer) = 0.      ! remove gas from the permeable layer
                   ENDIF
                   
                   IF ( i_imper_layer .EQ. 0 ) THEN ! ice-atmosphere bubble transfer
                      f_bub(jn) = f_bub(jn) + cbub_i_bio(jn,layer) *
     &                deltaz_i_bio(layer) / ddtb
                      cbub_i_bio(jn,layer) = 0.      ! remove gas from the permeable layer
                   ENDIF ! i_imper_layer
                 
                 END DO
                
                ENDIF ! flag_active, biotr_i_typ
                
             END DO ! jn
         
      ENDIF ! ln_bubrise

!
!------------------------------------------------------------------------------!
! 4) Update total gas concentration
!------------------------------------------------------------------------------!
!

      DO jn = 1, ntra_bio
         
         IF ( ( flag_active(jn)            ) .AND. 
     &        ( biotr_i_typ(jn) .EQ. 'gas' ) .AND.
     &        ln_write_gas ) THEN

         DO layer = 1, nlay_bio  ! dégazage des couches 2 à nlay_bio
            c_gtot_i(jn,layer) = cbu_i_bio(jn,layer) + 
     &                           cbub_i_bio(jn,layer)  
         END DO

         WRITE(numout,*) '    *** After ice_gas *** '

         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 )
         WRITE(numout,*) ' cbub_i_bio : ', ( cbub_i_bio(jn,layer), ! concentration of tracers in bubbles in bulk ice
     &                                   layer = 1, nlay_bio )
         WRITE(numout,*) ' c_gtot_i : ', ( c_gtot_i(jn,layer),   ! concentration of tracers in bubbles + bulk ice
     &                                   layer = 1, nlay_bio )
         WRITE(numout,*)
         
         ENDIF
         
      ENDDO

!
!------------------------------------------------------------------------------!
! 4) Conservation test and output
!------------------------------------------------------------------------------!
!

      CALL ice_gas_column( ntra_bio , c_gtot_i, deltaz_i_bio, .TRUE.,  
     &                     mt_i_gas_final )

      IF ( ln_write_gas ) THEN
         DO jn = 1, ntra_bio
            WRITE(numout,*) ' tracer         : ', biotr_i_nam(jn)
            WRITE(numout,*) ' mt_i_gas_init  : ', mt_i_gas_init(jn)
            WRITE(numout,*) ' mt_i_gas_final : ', mt_i_gas_final(jn)
         END DO
      ENDIF

      DO jn = 1, ntra_bio
         zfbub = f_bub(jn) * ddtb
         WRITE(numout,*) ' tracer : ', biotr_i_nam(jn)
         WRITE(numout,*) ' zfbub  : ', zfbub

      IF ( ABS( mt_i_gas_init(jn) - mt_i_gas_final(jn) - zfbub ) 
     &                                            .GT. 1.0e-10 )
     &   THEN
         WRITE(numout,*)
         WRITE(numout,*) ' ALERT!!! TOTAL GAS NOT CONSERVED '
         WRITE(numout,*) '      at time step : ', numit
         WRITE(numout,*) '         tracer    : ', jn
         WRITE(numout,*) ' mt_i_gas_init     : ', mt_i_gas_init(jn)
         WRITE(numout,*) ' mt_i_gas_final    : ', mt_i_gas_final(jn)
         WRITE(numout,*) ' diff    : ', mt_i_gas_init(jn) - 
     &                                  mt_i_gas_final(jn) - zfbub
         WRITE(numout,*) ' zfbub   : ', zfbub
         WRITE(numout,*)
      ENDIF
      ENDDO
      
!
!------------------------------------------------------------------------------!
! x) Update DIC and carbonate chemistry
!------------------------------------------------------------------------------!
!
      IF ( ln_carbon ) CALL ice_carb_chem

!------------------------------------------------------------------------------!
! X) End of the routine
!------------------------------------------------------------------------------!
!
      WRITE(numout,*)
      WRITE(numout,*) ' End of ice_gas '
      WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
     

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

      END