Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/bio.com
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/bio.com	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/bio.com	(revision 20)
@@ -66,5 +66,9 @@
      &   jn_car,      !:  ...
      &   jn_co2,      !:  ...
-     &   jn_cal
+     &   jn_cal,      !:  ...
+     &   jn_aon,      !:  ...
+     &   jn_eon,      !:  ...
+     &   jn_aop,      !:  ...           
+     &   jn_eop
 
       REAL(8), DIMENSION(ntra_bio_max,nlay_bio) ::
@@ -118,5 +122,9 @@
      &   fdiff,
      &   csat_gas,
-     &   sol_gas
+     &   sol_gas,
+     &   FDSI_AD,
+     &   FDIN_AD,
+     &   FDIP_AD
+
 
       COMMON /bioini/
@@ -137,5 +145,9 @@
      &   jn_car,      !:  ...
      &   jn_co2,      !:  ...
-     &   jn_cal
+     &   jn_cal,      !:  ...
+     &   jn_aon,      !:  ...
+     &   jn_eon,      !:  ...
+     &   jn_aop,      !:
+     &   jn_eop
 
       COMMON /biophyparams/
@@ -200,6 +212,8 @@
      &   ln_gasflux              , !: activate gas flux to atm
      &   ln_carbon               , !: activate carbon cycle
-     &   ln_ikaite                 !: activate CaCO3 
-     
+     &   ln_ikaite               , !: activate CaCO3 
+     &   ln_decoupNC             , !: activate N cycle disconnected from C cycle   
+     &   ln_decoupPC
+
       INTEGER(4) ::
      &   nn_bio_opt              , !: type of biological model (0=NP, 1=NPD Redfield)
@@ -258,5 +272,7 @@
      &   e_thr_gasflux           ,
      &   sursat_gas              ,
-     &   caco3_time               
+     &   caco3_time              ,
+     &   ln_decoupNC             ,
+     &   ln_decoupPC
 
       !-----------------------
@@ -329,5 +345,22 @@
      &   lim_tem(nlay_bio)       , !: temperature limitation
      &   lim_sal(nlay_bio)       , !: salt limitation
-     &   chlC_bio(nlay_bio)        !: interactive chl-a / C ratio
+     &   chlC_bio(nlay_bio)      , !: interactive chl-a / C ratio
+     &   syn_N(nlay_bio)         , !:
+     &   lys_N(nlay_bio)         , !:
+     &   rsp_N(nlay_bio)         , !:
+     &   rem_N(nlay_bio)         , !:
+     &   syn_P(nlay_bio)         , !:
+     &   lys_P(nlay_bio)         , !:
+     &   rsp_P(nlay_bio)         , !:
+     &   rem_P(nlay_bio)         , !:
+     &   N_C_alg(nlay_bio)       , !: N/C ratio in algae (if ln_decoupNC = TRUE)
+     &   N_C_det(nlay_bio)       , !: N/C ratio in detritic matter (if ln_decoupNC = TRUE)
+     &   lim_din_stock(nlay_bio) , !: limitation en stock DIN
+     &   lim_dip_stock(nlay_bio) , !:
+     &   lim_dsi_stock(nlay_bio) , !:
+     &   P_C_alg(nlay_bio)       , !: P/C ratio in algae (if ln_decoupPC = TRUE)
+     &   P_C_det(nlay_bio)       , !: P/C ratio in algae (if ln_decoupPC = TRUE)
+     &   N_P_alg(nlay_bio)       , !: N/P ratio in algae (if ln_decoupPC = TRUE)
+     &   N_P_det(nlay_bio)         !: N/P ratio in algae (if ln_decoupPC = TRUE)
 
       COMMON /biosources/
@@ -345,5 +378,23 @@
      &   lim_tem                 , !:
      &   lim_sal                 , !:
-     &   chlC_bio                  !: 
+     &   chlC_bio                , !: 
+     &   syn_N                   , !:
+     &   lys_N                   , !:
+     &   rsp_N                   , !:
+     &   rem_N                   , !:
+     &   syn_P                   , !:
+     &   lys_P                   , !:
+     &   rsp_P                   , !:
+     &   rem_P                   , !:
+     &   N_C_alg                 , !:
+     &   N_C_det                 , !:
+     &   lim_din_stock           , !:
+     &   lim_dip_stock           , !:
+     &   lim_dsi_stock           , !:
+     &   P_C_alg                 , !:
+     &   P_C_det                 , !:
+     &   N_P_alg                 , !:
+     &   N_P_det                   !:
+
 
       COMMON /biorad/             
@@ -374,5 +425,13 @@
      &   po4_c                   , !: N cell quota in diatoms
      &   oxy_c                   , !: Oxygen cell quota in diatoms
-     &   c_molar                   !: carbon molar mass
+     &   c_molar                 , !: carbon molar mass
+     &   ksyn_N                  , !: intensity of N flux compared with C flux
+     &   klys_N                  , !: same
+     &   krsp_N                  , !: same
+     &   krem_N                  , !: same
+     &   ksyn_P                  , !: intensity of N flux compared with C flux
+     &   klys_P                  , !: same
+     &   krsp_P                  , !: same
+     &   krem_P                    !: same
 
       COMMON /bioparams/
@@ -381,4 +440,7 @@
      &   krem_bio, krsp_bio, 
      &   rg_bio, rg_bac, chla_c, Estar, si_c, no3_c, po4_c, oxy_c,
-     &   lim_sal_wid, lim_sal_smax, c_molar
-
+     &   lim_sal_wid, lim_sal_smax, c_molar, ksyn_N, klys_N, 
+     &   krsp_N, krem_N, ksyn_P, klys_P, krsp_P, 
+     &   krem_P
+ 
+
Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_ini.f
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_ini.f	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_ini.f	(revision 20)
@@ -62,5 +62,5 @@
 
       IF ( c_bio_model .EQ. 'KRILL' ) THEN
-         ntra_bio = 21
+         ntra_bio = 23
       ENDIF
 !
@@ -117,4 +117,8 @@
       jn_arg = 12 ! 17
       jn_cal = 13 ! 21
+      jn_aon = 14
+      jn_eon = 15
+      jn_aop = 16
+      jn_eop = 17
 
       ! Assign Tracer names
@@ -132,4 +136,8 @@
       biotr_i_nam(jn_oxy) = 'Oxy'//zblank
       biotr_i_nam(jn_nit) = 'Nit'//zblank
+      biotr_i_nam(jn_aon) = 'AoN'//zblank
+      biotr_i_nam(jn_eon) = 'eoN'//zblank
+      biotr_i_nam(jn_aop) = 'AoP'//zblank
+      biotr_i_nam(jn_eop) = 'eoP'//zblank
 
 !
@@ -167,5 +175,5 @@
       READ(numtra,*)
 
-      DO jn_read = 1, 13
+      DO jn_read = 1, 17
 
          READ(numtra,*)  c_read_name(jn_read), zchar, i_dummy1, zdummy2,
@@ -248,5 +256,5 @@
       flag_adsorb(jn) = .false.
       flag_diff(jn)   = .false.
-      nn_remp(jn)     = 1
+      nn_remp(jn)     = 2
       biotr_i_typ(jn) = 'prc'//zblank
       biotr_i_uni(jn) = 'mmol_m3'//zblank
@@ -259,5 +267,49 @@
       flag_adsorb(jn) = .false.
       flag_diff(jn)   = .false.
-      nn_remp(jn)     = 1
+      nn_remp(jn)     = 2
+      biotr_i_typ(jn) = 'prc'//zblank
+      biotr_i_uni(jn) = 'mmol_m3'//zblank
+
+      !--------------------------------------
+      ! AoN : Ice algae organic N
+      !--------------------------------------
+      jn = jn_aon
+      flag_active(jn) = .true.
+      flag_adsorb(jn) = .false.
+      flag_diff(jn)   = .false.
+      nn_remp(jn)     = 2
+      biotr_i_typ(jn) = 'prc'//zblank
+      biotr_i_uni(jn) = 'mmol_m3'//zblank
+
+      !--------------------------------------
+      ! eoN : Detritus organic N
+      !--------------------------------------
+      jn = jn_eon
+      flag_active(jn) = .true.
+      flag_adsorb(jn) = .false.
+      flag_diff(jn)   = .false.
+      nn_remp(jn)     = 2
+      biotr_i_typ(jn) = 'prc'//zblank
+      biotr_i_uni(jn) = 'mmol_m3'//zblank
+
+      !--------------------------------------
+      ! AoP : Ice algae organic P
+      !--------------------------------------
+      jn = jn_aop
+      flag_active(jn) = .true.
+      flag_adsorb(jn) = .false.
+      flag_diff(jn)   = .false.
+      nn_remp(jn)     = 2
+      biotr_i_typ(jn) = 'prc'//zblank
+      biotr_i_uni(jn) = 'mmol_m3'//zblank
+
+      !--------------------------------------
+      ! eoP : Detritus organic P
+      !--------------------------------------
+      jn = jn_eop
+      flag_active(jn) = .true.
+      flag_adsorb(jn) = .false.
+      flag_diff(jn)   = .false.
+      nn_remp(jn)     = 2
       biotr_i_typ(jn) = 'prc'//zblank
       biotr_i_uni(jn) = 'mmol_m3'//zblank
Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_param.f
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_param.f	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_param.f	(revision 20)
@@ -93,4 +93,8 @@
       READ(numbio,*) ln_ikaite   ! activate ikaite
       READ(numbio,*)
+      READ(numbio,*) ln_decoupNC ! N cycle diconnected from C cycle
+      READ(numbio,*)
+      READ(numbio,*) ln_decoupPC ! P cycle diconnected from C cycle      
+      READ(numbio,*)
       READ(numbio,*)
       READ(numbio,*)
@@ -147,4 +151,20 @@
       READ(numbio,*)
       READ(numbio,*) frem_bio    ! fraction of diatom loss that is remineralized in Si
+      READ(numbio,*)
+      READ(numbio,*) ksyn_N      ! intensity of N flux VS C flux
+      READ(numbio,*)
+      READ(numbio,*) klys_N
+      READ(numbio,*)
+      READ(numbio,*) krsp_N
+      READ(numbio,*)
+      READ(numbio,*) krem_N
+      READ(numbio,*)
+      READ(numbio,*) ksyn_P
+      READ(numbio,*)
+      READ(numbio,*) klys_P
+      READ(numbio,*)
+      READ(numbio,*) krsp_P
+      READ(numbio,*)
+      READ(numbio,*) krem_P
       READ(numbio,*)
       READ(numbio,*)
@@ -219,4 +239,5 @@
       WRITE(numout,*) ' ln_carbon : ', ln_carbon
       WRITE(numout,*) ' ln_ikaite : ', ln_ikaite
+      WRITE(numbio,*) ' ln_decoupNC:', ln_decoupNC
       WRITE(numout,*)
       WRITE(numout,*) ' nn_init   : ', nn_init
@@ -234,4 +255,8 @@
       WRITE(numout,*) ' cco2_ini  : ', cco2_ini
       WRITE(numout,*) ' ccal_ini  : ', ccal_ini
+      WRITE(numbio,*) ' caon_ini  : ', caon_ini
+      WRITE(numbio,*) ' ceon_ini  : ', ceon_ini
+      WRITE(numbio,*) ' caop_ini  : ', caop_ini
+      WRITE(numbio,*) ' ceop_ini  : ', ceop_ini
       WRITE(numout,*) ' cno3_w    : ', cno3_w
       WRITE(numout,*) ' cpo4_w    : ', cpo4_w
@@ -247,4 +272,6 @@
       WRITE(numout,*) ' cco2_w    : ', cco2_w
       WRITE(numout,*) ' ccal_w    : ', ccal_w
+      WRITE(numbio,*) ' caon_w    : ', caon_w
+      WRITE(numbio,*) ' ceon_w    : ', ceon_w
       WRITE(numout,*) ' mr_oxy    : ', mr_oxy
       WRITE(numout,*) ' mr_nit    : ', mr_nit
@@ -271,5 +298,13 @@
       WRITE(numout,*) ' krsp_bio  : ', krsp_bio 
       WRITE(numout,*) ' krem_bio  : ', krem_bio 
-      WRITE(numout,*) ' frem_bio  : ', frem_bio 
+      WRITE(numout,*) ' frem_bio  : ', frem_bio
+      WRITE(numout,*) ' ksyn_N    : ', ksyn_N
+      WRITE(numout,*) ' klys_N    : ', klys_N
+      WRITE(numout,*) ' krsp_N    : ', krsp_N
+      WRITE(numout,*) ' krem_N    : ', krem_N
+      WRITE(numout,*) ' ksyn_P    : ', ksyn_P
+      WRITE(numout,*) ' klys_P    : ', klys_P
+      WRITE(numout,*) ' krsp_P    : ', krsp_P
+      WRITE(numout,*) ' krem_P    : ', krem_P 
       WRITE(numout,*) ' nn_phs    : ', nn_phs  
       WRITE(numout,*) ' ek_bio    : ', ek_bio
Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_sms.f
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_sms.f	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_bio_sms.f	(revision 20)
@@ -82,4 +82,16 @@
 !-----------------------------------------------------------------------------!
 !
+      consN = 0
+
+      DO layer = layer_00, nlay_bio 
+         
+         trucN = cbu_i_bio(jn_din,layer) + cbu_i_bio(jn_eon,layer)
+     &              + cbu_i_bio(jn_aon,layer) + consN
+         consN    = trucN        
+
+      END DO
+
+      WRITE(numout,*) ' Conservation BIO N :', consN
+
 
       ! Reference of the upper layer
@@ -246,36 +258,38 @@
       END DO
 
-      ! Remineralization 
-      IF ( ln_rem ) THEN
-
-         SELECT CASE ( nn_bio_opt ) 
-
-         CASE(0) ! NP
-
-            rem_bio(layer_00:nlay_bio) = frem_bio * 
-     &                                   lys_bio(layer_00:nlay_bio)
-
-         CASE(1) ! NPDr
-
-            IF ( nn_rem .GE. 1 )   ! no T-dependence
-     &         rem_bio(layer_00:nlay_bio) = krem_bio * 
-     &                          cbu_i_bio(jn_eoc,layer_00:nlay_bio)
-
-         
-            IF ( nn_rem .EQ. 2 )   ! T-dependence
-     &         rem_bio(layer_00:nlay_bio) = rem_bio(layer_00:nlay_bio) *
-     &                       EXP( rg_bac * tc_bio(layer_00:nlay_bio) )
-
-         END SELECT
-
-      ENDIF
-
+      DO layer = layer_00, nlay_bio
+
+         IF ( ln_rem ) THEN
+
+            SELECT CASE ( nn_bio_opt ) 
+
+               CASE(0) ! NP
+
+                  rem_bio(layer) = frem_bio * 
+     &                             lys_bio(layer)
+
+               CASE(1) ! NPD
+ 
+                  IF ( nn_rem .EQ. 1 )   ! no T-dependence
+     &               rem_bio(layer) = krem_bio * 
+     &                                cbu_i_bio(jn_eoc,layer)
+
+                  IF ( nn_rem .EQ. 2 )   ! T-dependence
+     &               rem_bio(layer) = krem_bio * cbu_i_bio(jn_eoc,layer)
+     &                                * EXP( rg_bac * tc_bio(layer) )
+     
+            END SELECT
+
+         ENDIF
+
+      END DO
 
       !------------------------------------------------------------------------
       ! 2.3 Diatom Synthesis
       !------------------------------------------------------------------------
+      zeps = 0.000000001
 
       DO layer = layer_00, nlay_bio
-   
+
          IF ( ln_syn ) THEN
             ! raw synthesis (in s-1)
@@ -287,15 +301,34 @@
             zsyn1 = z_syn * cbu_i_bio(jn_aoc,layer) ! normal carbon production rate (mmolC/m3/s)
             zsyn2 = zsyn1;    zsyn3 = zsyn1;    zsyn4 = zsyn1
+
    
             IF ( ln_lim_dsi )                  ! prevent exhaustion of Si
      &         zsyn2 = cbu_i_bio(jn_dsi,layer) / ddtb / si_c - zeps 
   
-            IF ( ln_lim_no3 )                  ! prevent exhaustion of N
-     &         zsyn3 = cbu_i_bio(jn_din,layer) / ddtb / no3_c - zeps
-  
-            IF ( ln_lim_po4 )                  ! prevent exhaustion of P
-     &         zsyn4 = cbu_i_bio(jn_dip,layer) / ddtb / po4_c - zeps
-   
-            ! nutrient stock-limited PP rate
+            IF (ln_decoupNC) THEN ! C aussi limité en stock par DIN
+               zC  = MAX (cbu_i_bio(jn_aoc,layer) , zeps )
+               zN  = MAX (cbu_i_bio(jn_aon,layer) , zeps )
+               N_C_alg(layer)  = zN / zC
+!              N_C_alg(layer) = 0.14 
+               zsyn3 = cbu_i_bio(jn_din,layer) / (ddtb*N_C_alg(layer))
+     &                 - zeps
+            ELSE 
+               zsyn3 = cbu_i_bio(jn_din,layer) / (ddtb * no3_c) - zeps
+            ENDIF
+
+
+            IF (ln_decoupPC) THEN ! C aussi limité en stock par DIP
+               zC  = MAX (cbu_i_bio(jn_aoc,layer) , zeps )
+               zP  = MAX (cbu_i_bio(jn_aop,layer) , zeps )
+               P_C_alg(layer)  = zP / zC
+!              P_C_alg(layer) = 0.0094 
+               zsyn4 = cbu_i_bio(jn_dip,layer)/(ddtb * P_C_alg(layer))
+     &                 - zeps
+            ELSE
+               zsyn4 = cbu_i_bio(jn_dip,layer) / ddtb / po4_c - zeps
+            ENDIF
+
+
+             ! nutrient stock-limited PP rate
             syn_bio(layer) = MIN( zsyn1, zsyn2, zsyn3, zsyn4 ) 
 
@@ -312,11 +345,126 @@
    
 
-      END DO ! layer
-
-      !------------------------------------------------------------------------
-      ! 2.4 Update reservoirs
-      !------------------------------------------------------------------------
-
-      DO layer = layer_00, nlay_bio
+         ! -------------------------
+         ! -- Limitations en stocks
+         ! ------------------------
+
+         IF (ln_decoupNC) THEN
+             zzdin = syn_bio(layer) / (((1/N_C_alg(layer)) *
+     &               (cbu_i_bio(jn_din,layer)/ddtb)) - zeps)
+
+         ELSE 
+             zzdin = syn_bio(layer) / (((1/no3_c) *
+     &               (cbu_i_bio(jn_din,layer)/ddtb)) - zeps)
+         ENDIF
+
+
+         IF (ln_decoupPC) THEN
+             zzpo4 = syn_bio(layer) / (((1/P_C_alg(layer)) *
+     &               (cbu_i_bio(jn_dip,layer)/ddtb)) - zeps)
+
+         ELSE
+             zzpo4 = syn_bio(layer) / (((1/po4_c) *
+     &               (cbu_i_bio(jn_dip,layer)/ddtb)) - zeps)
+         ENDIF
+
+
+
+         zzsi  = syn_bio(layer) / (((1/si_c) *
+     &           (cbu_i_bio(jn_dsi,layer)/ddtb)) - zeps)
+
+
+         IF ( zzdin .EQ. 1 ) THEN
+            lim_din_stock(layer) = 1   !DIN est limitant en stocks
+         ELSE
+            lim_din_stock(layer) = 0   !DIN PAS limitant en stocks
+         ENDIF
+
+         IF ( zzpo4 .EQ. 1 ) THEN
+            lim_dip_stock(layer) = 1
+         ELSE
+            lim_dip_stock(layer) = 0
+         ENDIF
+
+         IF ( zzsi .EQ. 1) THEN
+            lim_dsi_stock(layer) = 1
+         ELSE
+            lim_dsi_stock(layer) = 0
+         ENDIF
+
+
+      !------------------------------------------------------------------------
+      ! N cycle 
+      !------------------------------------------------------------------------
+
+
+         IF ( ln_decoupNC ) THEN
+
+            !--- Nitrogen carbon algal ratio
+            zC    = MAX( cbu_i_bio(jn_aoc,layer) , zeps )
+            zN    = MAX( cbu_i_bio(jn_aon,layer) , zeps )
+            N_C_alg(layer)  = zN / zC
+!           zN_C  = MIN( cbu_i_bio(jn_aon,layer) / zC , N_C_max )
+!           zN_C  = MAX( N_C_min , zN_C )
+!           N_C_alg(layer) = zN_C  ! => marche pas
+!            N_C_alg(layer) = 0.14
+
+            !--- Nitrogen carbon detritic ratio
+            zC    = MAX( cbu_i_bio(jn_eoc,layer), zeps )
+            zNd   = MAX( cbu_i_bio(jn_eon,layer), zeps ) 
+            N_C_det(layer) = zNd / zC
+!            N_C_det(layer) = 0.14
+
+            !--- Process rates
+            syn_N(layer) = N_C_alg(layer) * ksyn_N * syn_bio(layer)
+            rsp_N(layer) = N_C_alg(layer) * krsp_N * rsp_bio(layer)
+            lys_N(layer) = N_C_alg(layer) * klys_N * lys_bio(layer)
+            rem_N(layer) = N_C_det(layer) * krem_N * rem_bio(layer)
+
+
+         ENDIF ! ln_decoupNC
+
+      !------------------------------------------------------------------------
+      ! P cycle 
+      !------------------------------------------------------------------------
+
+
+         IF ( ln_decoupPC ) THEN
+
+            !--- P carbon algal ratio
+            zC    = MAX( cbu_i_bio(jn_aoc,layer) , zeps )
+            zP    = MAX( cbu_i_bio(jn_aop,layer) , zeps )
+            P_C_alg(layer)  = zP / zC
+
+            !--- P carbon detritic ratio
+            zC    = MAX( cbu_i_bio(jn_eoc,layer), zeps )
+            zPd   = MAX( cbu_i_bio(jn_eop,layer), zeps )
+            P_C_det(layer) = zPd / zC
+
+            !--- Process rates
+            syn_P(layer) = P_C_alg(layer) * ksyn_P * syn_bio(layer)
+            rsp_P(layer) = P_C_alg(layer) * krsp_P * rsp_bio(layer)
+            lys_P(layer) = P_C_alg(layer) * klys_P * lys_bio(layer)
+            rem_P(layer) = P_C_det(layer) * krem_P * rem_bio(layer)
+
+         ENDIF ! ln_decoupPC
+
+
+         zPa   = MAX( cbu_i_bio(jn_aop,layer) , zeps )
+         zNt   = MAX( cbu_i_bio(jn_aon,layer) , zeps )
+         N_P_alg(layer)  = zNt / zPa
+         zPe   = MAX( cbu_i_bio(jn_eop,layer) , zeps )
+         zNtd  = MAX( cbu_i_bio(jn_eon,layer) , zeps )
+         N_P_det(layer)  = zNtd / zPe
+
+
+
+      !------------------------------------------------------------------------
+      ! 2.4 Update reservoirs (including N cycle reservoirs)
+      !------------------------------------------------------------------------
+
+        ! -------------------------------------------------------
+        ! -- 1) Increments
+        ! ------------------------------------------------------
+
 
          ! DSi uptake
@@ -325,16 +473,4 @@
          zd_dsi2 = - cbu_i_bio(jn_dsi,layer)        ! Maximum DSi uptake
          zd_dsi  = MAX ( zd_dsi1, zd_dsi2 )
-
-         ! NO3 uptake
-         zd_no31 = no3_c * ( rem_bio(layer) +
-     &             rsp_bio(layer) - syn_bio(layer) ) * ddtb
-         zd_no32 = - cbu_i_bio(jn_din,layer)        ! Max NO3 uptake
-         zd_no3  = MAX ( zd_no31, zd_no32 )
-
-         ! update PO4
-         zd_po41 = po4_c * ( rem_bio(layer) +
-     &             rsp_bio(layer) - syn_bio(layer) ) * ddtb
-         zd_po42 = - cbu_i_bio(jn_dip,layer)        ! Max PO4 uptake
-         zd_po4  = MAX ( zd_po41, zd_po42 )
 
          ! update algae
@@ -350,13 +486,117 @@
             zd_eoc2 = - cbu_i_bio(jn_eoc,layer)     ! Max detrital carbon loss
             zd_eoc  = MAX( zd_eoc1, zd_eoc2 )
-         ENDIF
+         ENDIF ! nn_bio_opt
+
+
+         IF (ln_decoupNC) THEN
+
+            !--- Increments
+            zd_din1 = ( rsp_N(layer) - syn_N(layer) +
+     &                  rem_N(layer) ) * ddtb
+            zd_din2 = - cbu_i_bio(jn_din,layer) !PEU PAS EPUISER PLUS
+                                                !QUE CE QUIL Y A !!!!
+            zd_din  = MAX ( zd_din1, zd_din2 )  ! interessant dans le cas
+                                                ! ou zd_din1 est négatif
+                                                ! zd_din1 peut pas être
+                                                ! plus petit que zd_din2 
+
+            zd_aon1 = ( syn_N(layer) - rsp_N(layer) - lys_N(layer) ) *
+     &                ddtb
+            zd_aon2 = - cbu_i_bio(jn_aon,layer)
+            zd_aon  = MAX ( zd_aon1, zd_aon2 )
+
+            zd_eon1 = ( lys_N(layer) - rem_N(layer) ) * ddtb
+            zd_eon2 = - cbu_i_bio(jn_eon,layer)
+            zd_eon  = MAX ( zd_eon1, zd_eon2 )
+
+            !--- Correct rates to prevent leaks
+            IF ( zd_din .EQ. zd_din2 )
+     &         syn_N(layer) = MIN( syn_N(layer), -zd_din2/ddtb)
+
+            IF ( zd_aon .EQ. zd_aon2 )
+     &         lys_N(layer) = MIN( lys_N(layer), -zd_aon2/ddtb)
+
+            IF ( zd_eon .EQ. zd_eon2 )
+     &         rem_N(layer) = MIN( rem_N(layer), -zd_eon2/ddtb)
+
+
+         ELSE
+            zd_din1 = no3_c * ( rem_bio(layer) +
+     &                rsp_bio(layer) - syn_bio(layer) ) * ddtb
+            zd_din2 = - cbu_i_bio(jn_din,layer)        ! Max NO3 uptake
+            zd_din  = MAX ( zd_din1, zd_din2 )     
+             
+         ENDIF ! ln_decoupNC
+
+!!!!!!!!!!!!!!!!!!
+
+         IF (ln_decoupPC) THEN
+
+            !--- Increments
+            zd_dip1 = ( rsp_P(layer) - syn_P(layer) +
+     &                  rem_P(layer) ) * ddtb
+            zd_dip2 = - cbu_i_bio(jn_dip,layer) !PEU PAS EPUISER PLUS
+                                                !QUE CE QUIL Y A !!!!
+            zd_dip  = MAX ( zd_dip1, zd_dip2 )  ! interessant dans le
+                                                ! ou zd_din1 est négatif
+                                                ! zd_din1 peut pas être
+                                                ! plus petit que zd_din2 
+
+            zd_aop1 = ( syn_P(layer) - rsp_P(layer) - lys_P(layer) ) *
+     &                ddtb
+            zd_aop2 = - cbu_i_bio(jn_aop,layer)
+            zd_aop  = MAX ( zd_aop1, zd_aop2 )
+
+            zd_eop1 = ( lys_P(layer) - rem_P(layer) ) * ddtb
+            zd_eop2 = - cbu_i_bio(jn_eop,layer)
+            zd_eop  = MAX ( zd_eop1, zd_eop2 )
+
+            !--- Correct rates to prevent leaks
+            IF ( zd_dip .EQ. zd_dip2 )
+     &         syn_P(layer) = MIN( syn_P(layer), -zd_dip2/ddtb)
+
+            IF ( zd_aop .EQ. zd_aop2 )
+     &         lys_P(layer) = MIN( lys_P(layer), -zd_aop2/ddtb)
+
+            IF ( zd_eop .EQ. zd_eop2 )
+     &         rem_P(layer) = MIN( rem_P(layer), -zd_eop2/ddtb)
+
+
+         ELSE
+            zd_dip1 = po4_c * ( rem_bio(layer) +
+     &                rsp_bio(layer) - syn_bio(layer) ) * ddtb
+            zd_dip2 = - cbu_i_bio(jn_dip,layer)        ! Max DIP uptake
+            zd_dip  = MAX ( zd_dip1, zd_dip2 )
+
+         ENDIF ! ln_decoupPC
+
+
+        ! -------------------------------------------------------
+        ! -- 2) Add Increments
+        ! ------------------------------------------------------
+
+         IF (ln_decoupNC) THEN
+            cbu_i_bio(jn_din,layer) = cbu_i_bio(jn_din,layer) + zd_din
+            cbu_i_bio(jn_aon,layer) = cbu_i_bio(jn_aon,layer) + zd_aon
+            cbu_i_bio(jn_eon,layer) = cbu_i_bio(jn_eon,layer) + zd_eon
+         ELSE
+            cbu_i_bio(jn_din,layer) = cbu_i_bio(jn_din,layer) + zd_din
+         ENDIF
+
+         IF (ln_decoupPC) THEN
+            cbu_i_bio(jn_dip,layer) = cbu_i_bio(jn_dip,layer) + zd_dip
+            cbu_i_bio(jn_aop,layer) = cbu_i_bio(jn_aop,layer) + zd_aop
+            cbu_i_bio(jn_eop,layer) = cbu_i_bio(jn_eop,layer) + zd_eop
+         ELSE
+            cbu_i_bio(jn_dip,layer) = cbu_i_bio(jn_dip,layer) + zd_dip
+         ENDIF
+
 
          cbu_i_bio(jn_dsi,layer) = cbu_i_bio(jn_dsi,layer) + zd_dsi
-         cbu_i_bio(jn_din,layer) = cbu_i_bio(jn_din,layer) + zd_no3
-         cbu_i_bio(jn_dip,layer) = cbu_i_bio(jn_dip,layer) + zd_po4
          cbu_i_bio(jn_aoc,layer) = cbu_i_bio(jn_aoc,layer) + zd_daf
 
          IF ( nn_bio_opt .GE. 1 ) 
      &      cbu_i_bio(jn_eoc,layer) = cbu_i_bio(jn_eoc,layer) + zd_eoc
+
 
          IF ( ln_write_bio ) THEN
@@ -375,4 +615,7 @@
             WRITE(numout,*) ' lim_tem   : ', lim_tem(layer)
             WRITE(numout,*) ' lim_sal   : ', lim_sal(layer)
+            WRITE(numout,*) ' lim_din_stock : ', lim_din_stock(layer)
+            WRITE(numout,*) ' lim_dip_stock : ', lim_dip_stock(layer)
+            WRITE(numout,*) ' lim_dsi_stock : ', lim_dsi_stock(layer)
             WRITE(numout,*)
             WRITE(numout,*) ' tc_bio    : ', tc_bio(layer)  
@@ -381,19 +624,38 @@
             WRITE(numout,*)
             WRITE(numout,*) ' zd_dsi    : ', zd_dsi
-            WRITE(numout,*) ' zd_no3    : ', zd_no3
             WRITE(numout,*) ' zd_po4    : ', zd_po4
             WRITE(numout,*) ' zd_daf    : ', zd_daf
             WRITE(numout,*) ' zd_eoc    : ', zd_eoc
+            WRITE(numout,*) ' zd_din    : ', zd_din
             WRITE(numout,*)
             WRITE(numout,*) ' dSi       : ', cbu_i_bio(jn_dsi,layer)
-            WRITE(numout,*) ' dIN       : ', cbu_i_bio(jn_din,layer)
+            WRITE(numout,*) ' dIN       : ', cbu_i_bio(jn_din,layer) 
             WRITE(numout,*) ' dIP       : ', cbu_i_bio(jn_dip,layer)
             WRITE(numout,*) ' AoC       : ', cbu_i_bio(jn_aoc,layer)
             WRITE(numout,*) ' eoC       : ', cbu_i_bio(jn_eoc,layer)
-            WRITE(numout,*)
+            WRITE(numout,*) ' krem_N    : ', krem_N
+
+        
+            IF (ln_decoupNC) THEN 
+               WRITE(numout,*) ' syn_N    : ', syn_N(layer)
+               WRITE(numout,*) ' rsp_N    : ', rsp_N(layer)
+               WRITE(numout,*) ' lys_N    : ', lys_N(layer)
+               WRITE(numout,*) ' rem_N    : ', rem_N(layer)
+               WRITE(numout,*)
+               WRITE(numout,*) ' zd_aon    : ', zd_aon
+               WRITE(numout,*) ' zd_eon    : ', zd_eon
+               WRITE(numout,*)
+               WRITE(numout,*) ' AoN       : ', cbu_i_bio(jn_aon,layer)
+               WRITE(numout,*) ' eoN       : ', cbu_i_bio(jn_eon,layer)
+               WRITE(numout,*)
+               WRITE(numout,*) ' N_C_alg   : ', N_C_alg(layer)
+               WRITE(numout,*) ' N_C_det   : ', N_C_det(layer)
+            ENDIF ! ln_decoupNC
 
          ENDIF ! ln_write_bio
 
       END DO ! layer
+   
+
 
       !------------------------------------------------------------------------
@@ -408,9 +670,9 @@
      &                  .AND. flag_active(jn_alk) ) THEN
 
-            zd_dic1   = - zncp * ddtb
+            zd_dic1   = - zncp * ddtb                  !! Bizarre
             zd_dic2   = - cbu_i_bio(jn_dic,layer)
             zd_dic    = MAX( zd_dic1 , zd_dic2 )
 
-            zd_alk1   =   no3_c * zncp * ddtb
+            zd_alk1   =   no3_c * zncp * ddtb           !! Bizarre
             zd_alk2   = - cbu_i_bio(jn_alk,layer)
             zd_alk    = MAX( zd_alk1 , zd_alk2 )
Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_gas.f
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_gas.f	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_gas.f	(revision 20)
@@ -48,4 +48,5 @@
 !
 
+
       !--- Gas concentrations
       DO jn = 1, ntra_bio
@@ -59,6 +60,10 @@
          ENDIF
 
-         ! Carbonate chemistry
-         IF ( ln_carbon ) CALL ice_carb_chem
+      END DO
+
+      ! Carbonate chemistry
+      IF ( ln_carbon ) CALL ice_carb_chem
+
+      DO jn = 1, ntra_bio
 
          IF ( ( flag_active(jn)            ) .AND.
Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_output.f
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_output.f	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_output.f	(revision 20)
@@ -211,4 +211,21 @@
      .                   "kgNaCl/dm3","time","-","-","-" )
          WRITE(numout,*) ' oce_sal dimension created '
+
+         CALL CF_CREATE_VAR("divflux_AD","Sal flux divergence",
+     .                   "...","time","z_i_p","-","-" )
+         WRITE(numout,*) ' divflux_AD created '
+
+         CALL CF_CREATE_VAR("FDSI_AD","DSI flux for Advection",
+     .                   "mmol/m2/s","time","-","-","-" )
+         WRITE(numout,*) ' FDSI_AD created '
+
+         CALL CF_CREATE_VAR("FDIN_AD","DIN flux for Advection",
+     .                   "mmol/m2/s","time","-","-","-" )
+         WRITE(numout,*) ' FDIN_AD created '
+
+         CALL CF_CREATE_VAR("FDIP_AD","DIP flux for Advection",
+     .                   "mmol/m2/s","time","-","-","-" )
+         WRITE(numout,*) ' FDIP_AD created '
+
 
          ! Forcing
@@ -421,4 +438,36 @@
          WRITE(numout,*) ' ChlC_bio created '
 
+        ! N_C_alg
+         CALL CF_CREATE_VAR("N_C_alg", "N_C_alg",
+     &        "mmol.m-3","time","z_i_b","-","-")
+         WRITE(numout,*) ' N_C_alg created '
+
+        ! N_C_det
+         CALL CF_CREATE_VAR("N_C_det", "N_C_det",
+     &        "mmol.m-3","time","z_i_b","-","-")
+         WRITE(numout,*) ' N_C_det created '
+
+        ! P_C_alg
+         CALL CF_CREATE_VAR("P_C_alg", "P_C_alg",
+     &        "mmol.m-3","time","z_i_b","-","-")
+         WRITE(numout,*) ' P_C_alg created '
+
+        ! P_C_det
+         CALL CF_CREATE_VAR("P_C_det", "P_C_det",
+     &        "mmol.m-3","time","z_i_b","-","-")
+         WRITE(numout,*) ' P_C_det created '
+
+        ! N_P_alg
+         CALL CF_CREATE_VAR("N_P_alg", "N_P_alg",
+     &        "mmol.m-3","time","z_i_b","-","-")
+         WRITE(numout,*) ' N_P_alg created '
+
+        ! N_P_det
+         CALL CF_CREATE_VAR("N_P_det", "N_P_det",
+     &        "mmol.m-3","time","z_i_b","-","-")
+         WRITE(numout,*) ' N_P_det created '
+
+
+
          CALL CF_CREATE_VAR("Chla_phy", "Chlorophyll a on phy grid",
      &        "mmol.m-3","time","z_i_p","-","-" )
@@ -492,5 +541,17 @@
      .                   "-","time","z_i_b","-","-" )
          WRITE(numout,*) ' Sal limitation created '
-         
+
+         CALL CF_CREATE_VAR("lim_din_stock","DIN stock limitation",
+     .                   "-","time","z_i_b","-","-" )
+         WRITE(numout,*) ' DIN stock limitation created '
+
+         CALL CF_CREATE_VAR("lim_dip_stock","DIP stock limitation",
+     .                   "-","time","z_i_b","-","-" )
+         WRITE(numout,*) ' DIP stock limitation created '
+
+         CALL CF_CREATE_VAR("lim_dsi_stock","DSI stock limitation",
+     .                   "-","time","z_i_b","-","-" )
+         WRITE(numout,*) ' DSI stock limitation created '
+
          CALL CF_CREATE_VAR("divf_si","Flux divergence silicate",
      .                   "mmol Si m-3 s-1","time","z_i_b","-","-" )
@@ -593,4 +654,17 @@
      &               1, 1, 1, REAL(oce_sal) )
 
+
+      CALL CF_WRITE (filenc, 'FDSI_AD', numit-nstart+1,
+     &               1, 1, 1, REAL(FDSI_AD) )
+
+      CALL CF_WRITE (filenc, 'FDIN_AD', numit-nstart+1,
+     &               1, 1, 1, REAL(FDIN_AD) )
+
+      CALL CF_WRITE (filenc, 'FDIP_AD', numit-nstart+1,
+     &               1, 1, 1, REAL(FDIP_AD) )
+
+
+
+
       ! Forcing
       !---------
@@ -703,4 +777,13 @@
       CALL CF_WRITE (filenc, 't_s ', numit-nstart+1, nlay_s, 1, 1, 
      &               dummy_vector )
+
+      ! sal divergence flux AD
+      DO layer = 1, nlay_i
+         dummy_vector(layer) = REAL(divflux_AD(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'divflux_AD', numit-nstart+1,
+     &               nlay_i, 1, 1, dummy_vector)
+
+
 
       ! Radiative transfer
@@ -916,4 +999,47 @@
      &               dummy_vector )
 
+      ! N_C_alg            
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(N_C_alg(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'N_C_alg', numit-nstart+1, nlay_i, 1, 1,
+     &               dummy_vector )
+
+      ! N_C_det            
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(N_C_det(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'N_C_det', numit-nstart+1, nlay_i, 1, 1,
+     &               dummy_vector )
+
+      ! P_C_alg            
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(P_C_alg(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'P_C_alg', numit-nstart+1, nlay_i, 1, 1,
+     &               dummy_vector )
+
+      ! P_C_det            
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(P_C_det(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'P_C_det', numit-nstart+1, nlay_i, 1, 1,
+     &               dummy_vector )
+
+      ! N_P_alg            
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(N_P_alg(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'N_P_alg', numit-nstart+1, nlay_i, 1, 1,
+     &               dummy_vector )
+
+      ! N_P_det            
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(N_P_det(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'N_P_det', numit-nstart+1, nlay_i, 1, 1,
+     &               dummy_vector )
+
+
       ! Chla, biological grid
       DO layer = 1, nlay_bio
@@ -1035,4 +1161,25 @@
       CALL CF_WRITE (filenc, 'lim_sal', numit-nstart+1, nlay_bio, 1, 1, 
      &               dummy_vector )
+
+
+      ! Stock limitation
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(lim_din_stock(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'lim_din_stock', numit-nstart+1, 
+     &               nlay_bio, 1, 1, dummy_vector )
+
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(lim_dip_stock(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'lim_dip_stock', numit-nstart+1, 
+     &               nlay_bio, 1, 1, dummy_vector )
+
+      DO layer = 1, nlay_bio
+         dummy_vector(layer) = REAL(lim_dsi_stock(layer))
+      END DO
+      CALL CF_WRITE (filenc, 'lim_dsi_stock', numit-nstart+1, 
+     &               nlay_bio, 1, 1, dummy_vector )
+
 
       ! Flux divergence, silicate
Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_sal_adv.f
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_sal_adv.f	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_sal_adv.f	(revision 20)
@@ -43,5 +43,4 @@
      &   zb                 ,    !: gravity drainage
      &   zc                 ,    !: upward advective flow
-     &   w_adv_br           ,    !: brine velocity
      &   zRae               ,    !: effective Ra
      &   ze                 ,    !: downward advective flow
@@ -49,4 +48,5 @@
      &   zindtbis           ,    !:
      &   zdiagbis                !:
+
 
       REAL(8), DIMENSION(nlay_i,3) :: !: dummy factors for tracer equation
@@ -75,5 +75,5 @@
      &   zlevel             ,    !: height of the water column
      &   zthdiff                 !: thermal diffusivity
-
+ 
       REAL(8), DIMENSION(nlay_i+1) ::  
      &   z_sbr_int               !: brine salinity at layer interfaces
@@ -300,4 +300,6 @@
       rho_br_GN = 1020.
       Rc_GN     = 1.01
+      !alpha_GN = 1.0e-3
+      !Rc_GN    = 1.0
 
       w_adv_br(:) = 0.0
@@ -322,5 +324,5 @@
 !------------------------------------------------------------------------------|
 !
-
+  
       DO layer = 1, nlay_i
          za(layer) = w_adv_br(layer) * ddtb / deltaz_i_phy(layer)
@@ -330,4 +332,5 @@
       sn_i_b(1) = z_sbr_i(1) * ( e_i_b(1) + za(1) ) + 
      &            z_sbr_i(2) * ( - za(1) ) 
+
 
       ! inner layers
@@ -342,8 +345,8 @@
      &                   z_sbr_i(nlay_i)   * ( e_i_b(nlay_i) + 
      &                   za(nlay_i)/2. ) - za(nlay_i) * oce_sal
+ 
 
       IF ( ln_write ) THEN
          WRITE(numout,*) 
-         WRITE(numout,*)
          WRITE(numout,*) ' sn_i_b   : ', ( sn_i_b(layer) , 
      &                   layer = 1, nlay_i )
@@ -384,4 +387,6 @@
 
       ENDIF ! ln_sal
+
+
 !
 !------------------------------------------------------------------------------|
Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_th.f
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_th.f	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_th.f	(revision 20)
@@ -105,10 +105,16 @@
       IF ( c_gravdr .EQ. 'DI' )
      &   CALL ice_sal_diff(nlay_i,1,1)                          ! Salt transport (Eddy-diffusion-based)
-
-      IF ( c_gravdr .EQ. 'AD' )
+      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' ) .AND. ln_trdiff) THEN   ! Bio transport
+         IF ( c_gravdr .EQ. 'DI' )
+     &      CALL ice_bio_diff(1,1,nlay_i) 
+      ENDIF
+
+      IF ( ( c_bio_model .EQ. 'KRILL' ) .AND. ln_trdiff ) THEN
+         IF ( c_gravdr .EQ. 'AD' )
+     &      CALL ice_bio_adv(1,1,nlay_i) 
+      ENDIF
 
       IF ( c_bio_model .EQ. 'KRILL' )
Index: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/thermo.com
===================================================================
--- branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/thermo.com	(revision 10)
+++ branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/thermo.com	(revision 20)
@@ -102,5 +102,6 @@
       common/ salt / beta_sal, s_i_new, s_i_snic, e_skel, q_summer,
      &               diff_br(maxnlay), rayleigh(maxnlay), fsb, fsbp,
-     &               w_flood, w_flush
+     &               w_flood, w_flush, w_adv_br(maxnlay), 
+     &               divflux_AD(maxnlay)
 
       ! units is a flux, including radab
