Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 10574)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 10622)
@@ -11,4 +11,5 @@
    !! 'key_medusa'          : MEDUSA model
    !! 'key_roam'            : MEDUSA extras for carbonate cycle
+   !! 'key_fabm'            : ERSEM model coupled via FABM
    !! 'key_karaml'          : Kara mixed layer depth
    !!---------------------------------------------------------------------------
@@ -38,5 +39,6 @@
       & ln_kara,            &
 #endif   
-      & hmld,               & 
+      & hmld,               &
+      & hmld_tref,          &
       & hmlp,               &
       & hmlpt 
@@ -45,19 +47,25 @@
       & c_asmbal,           &
       & nitbkg_r,           &
+      & nitavgbkg_r,        &
       & nitdin_r,           &
       & nitiaustr_r,        &
       & nitiaufin_r
+   USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
+      & asm_phyto2d_bal_medusa
+   USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC
+      & asm_phyto2d_bal_hadocc
+   USE asmphyto2dbal_ersem,  ONLY: & ! phyto2d balancing for ERSEM
+      & asm_phyto2d_bal_ersem
+   USE asmpco2bal, ONLY:    & ! pCO2 balancing
+      & asm_pco2_bal
 #if defined key_top
    USE trc, ONLY:           & ! passive tracer variables
       & trn,                &
-      & trb
+      & trb,                &
+      & nittrc000
    USE par_trc, ONLY:       & ! passive tracer parameters
       & jptra
 #endif
 #if defined key_medusa
-   USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
-      & asm_phyto2d_bal_medusa
-   USE asmpco2bal, ONLY:    & ! pCO2 balancing for MEDUSA
-      & asm_pco2_bal
    USE sms_medusa             ! MEDUSA parameters
    USE par_medusa             ! MEDUSA parameters
@@ -72,8 +80,4 @@
       & mld_max
 #elif defined key_hadocc
-   USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC
-      & asm_phyto2d_bal_hadocc
-   USE asmpco2bal, ONLY:    & ! pCO2 balancing for HadOCC
-      & asm_pco2_bal
    USE par_hadocc             ! HadOCC parameters
    USE had_bgc_const          ! HadOCC parameters
@@ -86,4 +90,11 @@
    USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio
       & cchl_p
+#elif defined key_fabm
+   USE par_fabm               ! FABM-ERSEM parameters
+   USE trc, ONLY:           & ! FABM-ERSEM diagnostic variables
+      & pgrow_avg,          &
+      & ploss_avg,          &
+      & phyt_avg,           &
+      & mld_max
 #endif
 
@@ -97,4 +108,6 @@
    PUBLIC  asm_bgc_init_bkg       ! called by asm_inc_init in asminc.F90
    PUBLIC  asm_bgc_bal_wri        ! called by nemo_gcm in nemogcm.F90
+   PUBLIC  asm_bgc_bkg_alloc      ! called by asm_bkg_wri in asmbkg.F90
+   PUBLIC  asm_bgc_bkg_tavg       ! called by asm_bkg_wri in asmbkg.F90
    PUBLIC  asm_bgc_bkg_wri        ! called by asm_bkg_wri in asmbkg.F90
    PUBLIC  phyto2d_asm_inc        ! called by bgc_asm_inc in asminc.F90
@@ -109,4 +122,7 @@
    LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom     log10(chlorophyll) increment
    LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_slchlnaninc = .FALSE. !: No surface nanophyto  log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_slchlpicinc = .FALSE. !: No surface picophyto  log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_slchldininc = .FALSE. !: No surface dinoflag   log10(chlorophyll) increment
    LOGICAL, PUBLIC :: ln_schltotinc  = .FALSE. !: No surface total      chlorophyll        increment
    LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total      log10(phyto C)     increment
@@ -146,4 +162,7 @@
    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchldia_bkginc ! slchldia inc
    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlnan_bkginc ! slchlnan inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlpic_bkginc ! slchlpic inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchldin_bkginc ! slchldin inc
    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: schltot_bkginc  ! schltot inc
    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphytot_bkginc ! slphytot inc
@@ -172,4 +191,18 @@
    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: chl_bkg         ! Background chl
    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: cchl_p_bkg      ! Background c:chl
+#if defined key_fabm
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: totalk_bkg      ! Background total alkalinity
+#endif
+
+   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: pgrow_avg_tavg  ! Time average of pgrow_avg
+   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: ploss_avg_tavg  ! Time average of ploss_avg
+   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: phyt_avg_tavg   ! Time average of phyt_avg
+   REAL(wp), DIMENSION(:,:,:,:), SAVE, ALLOCATABLE :: trn_tavg        ! Time average of trn
+#if defined key_hadocc
+   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: HADOCC_CHL_tavg ! Time average of HADOCC_CHL
+   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: cchl_p_tavg     ! Time average of cchl_p
+#elif defined key_fabm
+   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: totalk_tavg     ! Time average of O3_TA
+#endif
 
 #  include "domzgr_substitute.h90"
@@ -190,5 +223,5 @@
       !!------------------------------------------------------------------------
 
-#if ! defined key_top || ( ! defined key_hadocc && ! defined key_medusa )
+#if ! defined key_top || ( ! defined key_hadocc && ! defined key_medusa && ! defined key_fabm )
       CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', &
          &           ' but no compatible biogeochemical model is available' )
@@ -208,9 +241,16 @@
 #endif
 
+#if defined key_fabm
+      IF ( ln_pphinc ) THEN
+         CALL ctl_stop( ' Cannot currently assimilate pH into FABM-ERSEM' )
+      ENDIF
+#endif
+
       IF ( ( ln_phytobal ).AND.                                       &
          & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
-         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. &
-         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
-         & ( .NOT. ln_slphynoninc ) ) THEN
+         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. &
+         & ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. &
+         & ( .NOT. ln_schltotinc  ).AND.( .NOT. ln_slphytotinc ).AND. &
+         & ( .NOT. ln_slphydiainc ).AND.( .NOT. ln_slphynoninc ) ) THEN
          CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
             &           ' if not assimilating ocean colour,',                   &
@@ -221,5 +261,7 @@
       IF ( ( ln_balwri ).AND.                                         &
          & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
-         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. &
+         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. &
+         & ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. &
+         & ( .NOT. ln_schltotinc  ).AND.                              &
          & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
          & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
@@ -244,5 +286,6 @@
 
       IF ( ( ln_slchltotinc .OR. ln_schltotinc  ) .AND. &
-         & ( ln_slchldiainc .OR. ln_slchlnoninc ) ) THEN
+         & ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slchlnaninc .OR. &
+         &   ln_slchlpicinc .OR. ln_slchldininc ) ) THEN
          CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' )
       ENDIF
@@ -250,4 +293,9 @@
       IF ( ln_phytobal .AND.                                      &
          & ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. &
+#if defined key_fabm
+         &   (  ( ln_slchldiainc .OR. ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc ) .AND. &
+         &      ( ( .NOT. ln_slchldiainc ) .OR. ( .NOT. ln_slchlnaninc ) .OR. &
+         &        ( .NOT. ln_slchlpicinc ) .OR. ( .NOT. ln_slchldininc ) ) ) .OR. &
+#endif
          &   ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN
          CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
@@ -306,4 +354,19 @@
       ENDIF
       
+      IF ( ln_slchlnaninc ) THEN
+         ALLOCATE( slchlnan_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnan', slchlnan_bkginc )
+      ENDIF
+      
+      IF ( ln_slchlpicinc ) THEN
+         ALLOCATE( slchlpic_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlpic', slchlpic_bkginc )
+      ENDIF
+      
+      IF ( ln_slchldininc ) THEN
+         ALLOCATE( slchldin_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslchldin', slchldin_bkginc )
+      ENDIF
+      
       IF ( ln_schltotinc ) THEN
          ALLOCATE( schltot_bkginc(jpi,jpj) )
@@ -379,4 +442,5 @@
       
       IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
          & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
          & ln_slphynoninc ) THEN
@@ -516,6 +580,7 @@
       !!------------------------------------------------------------------------
 
-#if defined key_top && ( defined key_hadocc || defined key_medusa )
+#if defined key_top && ( defined key_hadocc || defined key_medusa || defined key_fabm )
       IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
          & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
          & ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN
@@ -557,4 +622,22 @@
          CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
          CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
+#elif defined key_fabm
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl1', tracer_bkg(:,:,:,jp_fabm_chl1) )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl2', tracer_bkg(:,:,:,jp_fabm_chl2) )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl3', tracer_bkg(:,:,:,jp_fabm_chl3) )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl4', tracer_bkg(:,:,:,jp_fabm_chl4) )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1c',  tracer_bkg(:,:,:,jp_fabm_p1c)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1n',  tracer_bkg(:,:,:,jp_fabm_p1n)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1p',  tracer_bkg(:,:,:,jp_fabm_p1p)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1s',  tracer_bkg(:,:,:,jp_fabm_p1s)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2c',  tracer_bkg(:,:,:,jp_fabm_p2c)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2n',  tracer_bkg(:,:,:,jp_fabm_p2n)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2p',  tracer_bkg(:,:,:,jp_fabm_p2p)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3c',  tracer_bkg(:,:,:,jp_fabm_p3c)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3n',  tracer_bkg(:,:,:,jp_fabm_p3n)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3p',  tracer_bkg(:,:,:,jp_fabm_p3p)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4c',  tracer_bkg(:,:,:,jp_fabm_p4c)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4n',  tracer_bkg(:,:,:,jp_fabm_p4n)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4p',  tracer_bkg(:,:,:,jp_fabm_p4p)  )
 #endif
          
@@ -588,4 +671,21 @@
             CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
             CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
+#elif defined key_fabm
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_z4c',  tracer_bkg(:,:,:,jp_fabm_z4c)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5c',  tracer_bkg(:,:,:,jp_fabm_z5c)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5n',  tracer_bkg(:,:,:,jp_fabm_z5n)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5p',  tracer_bkg(:,:,:,jp_fabm_z5p)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6c',  tracer_bkg(:,:,:,jp_fabm_z6c)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6n',  tracer_bkg(:,:,:,jp_fabm_z6n)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6p',  tracer_bkg(:,:,:,jp_fabm_z6p)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_n1p',  tracer_bkg(:,:,:,jp_fabm_n1p)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_n3n',  tracer_bkg(:,:,:,jp_fabm_n3n)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_n4n',  tracer_bkg(:,:,:,jp_fabm_n4n)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_n5s',  tracer_bkg(:,:,:,jp_fabm_n5s)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_o2o',  tracer_bkg(:,:,:,jp_fabm_o2o)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_o3c)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              )
+            totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
 #endif
          ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
@@ -596,4 +696,9 @@
             CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
             CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+#elif defined key_fabm
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_o3c)  )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
+            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              )
+            totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
 #endif
             CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
@@ -622,4 +727,9 @@
          CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
          CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+#elif defined key_fabm
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_o3c)  )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
+         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              )
+         totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
 #endif
          CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
@@ -687,4 +797,5 @@
 
          IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+            & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
             & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
             & ln_slphynoninc ) THEN
@@ -716,4 +827,38 @@
                CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
             ENDIF
+#elif defined key_fabm
+            CALL iom_rstput( kt, kt, inum, 'phy2d_chl1', phyto2d_balinc(:,:,:,jp_fabm_chl1) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_chl2', phyto2d_balinc(:,:,:,jp_fabm_chl2) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_chl3', phyto2d_balinc(:,:,:,jp_fabm_chl3) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_chl4', phyto2d_balinc(:,:,:,jp_fabm_chl4) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p1c',  phyto2d_balinc(:,:,:,jp_fabm_p1c)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p1n',  phyto2d_balinc(:,:,:,jp_fabm_p1n)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p1p',  phyto2d_balinc(:,:,:,jp_fabm_p1p)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p1s',  phyto2d_balinc(:,:,:,jp_fabm_p1s)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p2c',  phyto2d_balinc(:,:,:,jp_fabm_p2c)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p2n',  phyto2d_balinc(:,:,:,jp_fabm_p2n)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p2p',  phyto2d_balinc(:,:,:,jp_fabm_p2p)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p3c',  phyto2d_balinc(:,:,:,jp_fabm_p3c)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p3n',  phyto2d_balinc(:,:,:,jp_fabm_p3n)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p3p',  phyto2d_balinc(:,:,:,jp_fabm_p3p)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p4c',  phyto2d_balinc(:,:,:,jp_fabm_p4c)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p4n',  phyto2d_balinc(:,:,:,jp_fabm_p4n)  )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_p4p',  phyto2d_balinc(:,:,:,jp_fabm_p4p)  )
+            IF ( ln_phytobal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'phy2d_z4c',   phyto2d_balinc(:,:,:,jp_fabm_z4c)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_z5c',   phyto2d_balinc(:,:,:,jp_fabm_z5c)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_z5n',   phyto2d_balinc(:,:,:,jp_fabm_z5n)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_z5p',   phyto2d_balinc(:,:,:,jp_fabm_z5p)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_z6c',   phyto2d_balinc(:,:,:,jp_fabm_z6c)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_z6n',   phyto2d_balinc(:,:,:,jp_fabm_z6n)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_z6p',   phyto2d_balinc(:,:,:,jp_fabm_z6p)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_n1p',   phyto2d_balinc(:,:,:,jp_fabm_n1p)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_n3n',   phyto2d_balinc(:,:,:,jp_fabm_n3n)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_n4n',   phyto2d_balinc(:,:,:,jp_fabm_n4n)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_n5s',   phyto2d_balinc(:,:,:,jp_fabm_n5s)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_o2o',   phyto2d_balinc(:,:,:,jp_fabm_o2o)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_o3c',   phyto2d_balinc(:,:,:,jp_fabm_o3c)  )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_o3ba',  phyto2d_balinc(:,:,:,jp_fabm_o3ba) )
+            ENDIF
 #endif
          ENDIF
@@ -728,4 +873,22 @@
 #elif defined key_hadocc
             CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
+#elif defined key_fabm
+            CALL iom_rstput( kt, kt, inum, 'phy3d_chl1', phyto3d_balinc(:,:,:,jp_fabm_chl1) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_chl2', phyto3d_balinc(:,:,:,jp_fabm_chl2) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_chl3', phyto3d_balinc(:,:,:,jp_fabm_chl3) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_chl4', phyto3d_balinc(:,:,:,jp_fabm_chl4) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p1c',  phyto3d_balinc(:,:,:,jp_fabm_p1c)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p1n',  phyto3d_balinc(:,:,:,jp_fabm_p1n)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p1p',  phyto3d_balinc(:,:,:,jp_fabm_p1p)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p1s',  phyto3d_balinc(:,:,:,jp_fabm_p1s)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p2c',  phyto3d_balinc(:,:,:,jp_fabm_p2c)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p2n',  phyto3d_balinc(:,:,:,jp_fabm_p2n)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p2p',  phyto3d_balinc(:,:,:,jp_fabm_p2p)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p3c',  phyto3d_balinc(:,:,:,jp_fabm_p3c)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p3n',  phyto3d_balinc(:,:,:,jp_fabm_p3n)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p3p',  phyto3d_balinc(:,:,:,jp_fabm_p3p)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p4c',  phyto3d_balinc(:,:,:,jp_fabm_p4c)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p4n',  phyto3d_balinc(:,:,:,jp_fabm_p4n)  )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_p4p',  phyto3d_balinc(:,:,:,jp_fabm_p4p)  )
 #endif
          ENDIF
@@ -738,4 +901,7 @@
             CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
             CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#elif defined key_fabm
+            CALL iom_rstput( kt, kt, inum, 'pco2_o3c',  pco2_balinc(:,:,:,jp_fabm_o3c)  )
+            CALL iom_rstput( kt, kt, inum, 'pco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) )
 #endif
          ELSE IF ( ln_sfco2inc ) THEN
@@ -746,4 +912,7 @@
             CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
             CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#elif defined key_fabm
+            CALL iom_rstput( kt, kt, inum, 'fco2_o3c',  pco2_balinc(:,:,:,jp_fabm_o3c)  )
+            CALL iom_rstput( kt, kt, inum, 'fco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) )
 #endif
          ENDIF
@@ -756,4 +925,7 @@
             CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) )
             CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) )
+#elif defined key_fabm
+            CALL iom_rstput( kt, kt, inum, 'ph_o3c',  ph_balinc(:,:,:,jp_fabm_o3c)  )
+            CALL iom_rstput( kt, kt, inum, 'ph_o3ba', ph_balinc(:,:,:,jp_fabm_o3ba) )
 #endif
          ENDIF
@@ -774,5 +946,73 @@
    !!===========================================================================
 
-   SUBROUTINE asm_bgc_bkg_wri( kt, knum )
+   SUBROUTINE asm_bgc_bkg_alloc
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_bkg_alloc  ***
+      !!
+      !! ** Purpose :   allocate time-average arrays for background
+      !!
+      !! ** Method  :   allocate time-average arrays for background
+      !!
+      !! ** Action  :   allocate time-average arrays for background
+      !!
+      !! References :   asm_bkg_wri
+      !!------------------------------------------------------------------------
+      !!
+      INTEGER :: ierror
+      !!
+      !!------------------------------------------------------------------------
+      
+      ALLOCATE( pgrow_avg_tavg(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate pgrow_avg_tavg' )
+      ENDIF
+      pgrow_avg_tavg(:,:) = 0.0
+      
+      ALLOCATE( ploss_avg_tavg(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate ploss_avg_tavg' )
+      ENDIF
+      ploss_avg_tavg(:,:) = 0.0
+      
+      ALLOCATE( phyt_avg_tavg(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate phyt_avg_tavg' )
+      ENDIF
+      phyt_avg_tavg(:,:) = 0.0
+      
+      ALLOCATE( trn_tavg(jpi,jpj,jpk,jptra), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate trn_tavg' )
+      ENDIF
+      trn_tavg(:,:,:,:) = 0.0
+
+#if defined key_hadocc
+      ALLOCATE( HADOCC_CHL_tavg(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate HADOCC_CHL_tavg' )
+      ENDIF
+      HADOCC_CHL_tavg(:,:,:) = 0.0
+
+      ALLOCATE( cchl_p_tavg(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate cchl_p_tavg' )
+      ENDIF
+      cchl_p_tavg(:,:,:) = 0.0
+
+#elif defined key_fabm
+      ALLOCATE( totalk_tavg(jpi,jpj,jpk), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate totalk_tavg' )
+      ENDIF
+      totalk_tavg(:,:,:) = 0.0
+#endif
+                                 
+   END SUBROUTINE asm_bgc_bkg_alloc
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE asm_bgc_bkg_tavg(kt, pnumtimes_tavg)
       !!------------------------------------------------------------------------
       !!                    ***  ROUTINE asm_bgc_bkg_wri  ***
@@ -787,42 +1027,150 @@
       !!------------------------------------------------------------------------
       !!
-      INTEGER, INTENT(in   ) :: kt     ! Current time-step
-      INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file
-      !!
-      !!------------------------------------------------------------------------
-
+      INTEGER,  INTENT(in   ) :: kt               ! Current time-step
+      !!
+      REAL(wp), INTENT(in   ) :: pnumtimes_tavg   ! No of times to average over
+      !!
+      !!------------------------------------------------------------------------
+      
+      IF (kt == nittrc000) THEN
+         pgrow_avg_tavg(:,:) = 0.0
+         ploss_avg_tavg(:,:) = 0.0
+         phyt_avg_tavg(:,:)  = 0.0
+         trn_tavg(:,:,:,:)   = 0.0
 #if defined key_hadocc
-      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg             )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg             )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg              )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max               )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo'  , trn(:,:,:,jp_had_zoo) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn'  , trn(:,:,:,jp_had_pdn) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic'  , trn(:,:,:,jp_had_dic) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk'  , trn(:,:,:,jp_had_alk) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl'  , HADOCC_CHL(:,:,:)     )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         )
+         HADOCC_CHL_tavg(:,:,:) = 0.0
+         cchl_p_tavg(:,:,:)     = 0.0
+#elif defined key_fabm
+         totalk_tavg(:,:,:)  = 0.0
+#endif
+      ENDIF
+      
+      pgrow_avg_tavg(:,:) = pgrow_avg_tavg(:,:) + pgrow_avg(:,:) / pnumtimes_tavg
+      ploss_avg_tavg(:,:) = ploss_avg_tavg(:,:) + ploss_avg(:,:) / pnumtimes_tavg
+      phyt_avg_tavg(:,:)  = phyt_avg_tavg(:,:)  + phyt_avg(:,:)  / pnumtimes_tavg
+      trn_tavg(:,:,:,:)   = trn_tavg(:,:,:,:)   + trn(:,:,:,:)   / pnumtimes_tavg
+#if defined key_hadocc
+      HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL_tavg(:,:,:) + HADOCC_CHL(:,:,:) / pnumtimes_tavg
+      cchl_p_tavg(:,:,:)     = cchl_p_tavg(:,:,:)     + cchl_p(:,:,:)     / pnumtimes_tavg
+#elif defined key_fabm
+      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) + &
+         &                  fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg
+      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:)
+#endif
+                                 
+   END SUBROUTINE asm_bgc_bkg_tavg
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE asm_bgc_bkg_wri( kt, knum, ld_avgbkg )
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_bkg_wri  ***
+      !!
+      !! ** Purpose :   write out bgc background
+      !!
+      !! ** Method  :   write out bgc background
+      !!
+      !! ** Action  :   write out bgc background
+      !!
+      !! References :   asm_bkg_wri
+      !!------------------------------------------------------------------------
+      !!
+      INTEGER, INTENT(in   ) :: kt          ! Current time-step
+      INTEGER, INTENT(in   ) :: knum        ! i/o unit of increments file
+      LOGICAL, INTENT(in   ) :: ld_avgbkg   ! Averaged background?
+      !!
+      INTEGER                :: nitbgcbkg_r ! Period referenced to nit000
+      !!
+      !!------------------------------------------------------------------------
+      
+      IF (ld_avgbkg) THEN
+         nitbgcbkg_r = nitavgbkg_r
+      ELSE
+         nitbgcbkg_r = nitbkg_r
+         pgrow_avg_tavg(:,:) = pgrow_avg(:,:)
+         ploss_avg_tavg(:,:) = ploss_avg(:,:)
+         phyt_avg_tavg(:,:)  = phyt_avg(:,:)
+         trn_tavg(:,:,:,:)   = trn(:,:,:,:)
+#if defined key_hadocc
+         HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL(:,:,:)
+         cchl_p_tavg(:,:,:)     = cchl_p(:,:,:)
+#elif defined key_fabm
+         totalk_tavg(:,:,:)  = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
+         totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:)
+#endif
+      ENDIF
+
+#if defined key_hadocc
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg             )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg             )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg              )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max                    )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_nut'  , trn_tavg(:,:,:,jp_had_nut) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_phy'  , trn_tavg(:,:,:,jp_had_phy) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_zoo'  , trn_tavg(:,:,:,jp_had_zoo) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_pdn'  , trn_tavg(:,:,:,jp_had_pdn) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_dic'  , trn_tavg(:,:,:,jp_had_dic) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_alk'  , trn_tavg(:,:,:,jp_had_alk) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_chl'  , HADOCC_CHL(:,:,:)     )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         )
 #elif defined key_medusa
-      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg        )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg        )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg         )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max          )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd'  , trn(:,:,:,jpchd) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn'  , trn(:,:,:,jpphn) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd'  , trn(:,:,:,jpphd) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds'  , trn(:,:,:,jppds) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi'  , trn(:,:,:,jpzmi) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme'  , trn(:,:,:,jpzme) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din'  , trn(:,:,:,jpdin) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil'  , trn(:,:,:,jpsil) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer'  , trn(:,:,:,jpfer) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det'  , trn(:,:,:,jpdet) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc'  , trn(:,:,:,jpdtc) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic'  , trn(:,:,:,jpdic) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk'  , trn(:,:,:,jpalk) )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy'  , trn(:,:,:,jpoxy) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg        )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg        )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg         )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max               )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chn'  , trn_tavg(:,:,:,jpchn) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chd'  , trn_tavg(:,:,:,jpchd) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phn'  , trn_tavg(:,:,:,jpphn) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phd'  , trn_tavg(:,:,:,jpphd) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_pds'  , trn_tavg(:,:,:,jppds) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zmi'  , trn_tavg(:,:,:,jpzmi) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zme'  , trn_tavg(:,:,:,jpzme) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_din'  , trn_tavg(:,:,:,jpdin) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_sil'  , trn_tavg(:,:,:,jpsil) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_fer'  , trn_tavg(:,:,:,jpfer) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_det'  , trn_tavg(:,:,:,jpdet) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dtc'  , trn_tavg(:,:,:,jpdtc) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dic'  , trn_tavg(:,:,:,jpdic) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_alk'  , trn_tavg(:,:,:,jpalk) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_oxy'  , trn_tavg(:,:,:,jpoxy) )
+#elif defined key_fabm
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg               )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg               )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg                )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max                      )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl1'  , trn_tavg(:,:,:,jp_fabm_chl1) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl2'  , trn_tavg(:,:,:,jp_fabm_chl2) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl3'  , trn_tavg(:,:,:,jp_fabm_chl3) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl4'  , trn_tavg(:,:,:,jp_fabm_chl4) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1c'   , trn_tavg(:,:,:,jp_fabm_p1c)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1n'   , trn_tavg(:,:,:,jp_fabm_p1n)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1p'   , trn_tavg(:,:,:,jp_fabm_p1p)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1s'   , trn_tavg(:,:,:,jp_fabm_p1s)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2c'   , trn_tavg(:,:,:,jp_fabm_p2c)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2n'   , trn_tavg(:,:,:,jp_fabm_p2n)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2p'   , trn_tavg(:,:,:,jp_fabm_p2p)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3c'   , trn_tavg(:,:,:,jp_fabm_p3c)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3n'   , trn_tavg(:,:,:,jp_fabm_p3n)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3p'   , trn_tavg(:,:,:,jp_fabm_p3p)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4c'   , trn_tavg(:,:,:,jp_fabm_p4c)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4n'   , trn_tavg(:,:,:,jp_fabm_p4n)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4p'   , trn_tavg(:,:,:,jp_fabm_p4p)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z4c'   , trn_tavg(:,:,:,jp_fabm_z4c)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5c'   , trn_tavg(:,:,:,jp_fabm_z5c)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5n'   , trn_tavg(:,:,:,jp_fabm_z5n)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5p'   , trn_tavg(:,:,:,jp_fabm_z5p)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6c'   , trn_tavg(:,:,:,jp_fabm_z6c)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6n'   , trn_tavg(:,:,:,jp_fabm_z6n)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6p'   , trn_tavg(:,:,:,jp_fabm_z6p)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n1p'   , trn_tavg(:,:,:,jp_fabm_n1p)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n3n'   , trn_tavg(:,:,:,jp_fabm_n3n)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n4n'   , trn_tavg(:,:,:,jp_fabm_n4n)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n5s'   , trn_tavg(:,:,:,jp_fabm_n5s)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o2o'   , trn_tavg(:,:,:,jp_fabm_o2o)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3c'   , trn_tavg(:,:,:,jp_fabm_o3c)  )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ba'  , trn_tavg(:,:,:,jp_fabm_o3ba) )
+      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ta'  , totalk_tavg                  )
 #endif
 
@@ -911,4 +1259,13 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs
       REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
+#elif defined key_fabm
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldin ! Local chldin incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldin ! Local chldin bkg
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnan ! Local chlnan incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnan ! Local chlnan bkg
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlpic ! Local chlpic incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlpic ! Local chlpic bkg
 #endif
       !!------------------------------------------------------------------------
@@ -926,4 +1283,9 @@
 #elif defined key_hadocc
             zbkg_chltot(:,:) = chl_bkg(:,:,1)
+#elif defined key_fabm
+            zbkg_chltot(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1) + &
+               &               tracer_bkg(:,:,1,jp_fabm_chl2) + &
+               &               tracer_bkg(:,:,1,jp_fabm_chl3) + &
+               &               tracer_bkg(:,:,1,jp_fabm_chl4)
 #endif
             CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot )
@@ -934,8 +1296,12 @@
          ENDIF
 
-#if defined key_medusa
+#if defined key_medusa || defined key_fabm
          ! Diatom chlorophyll
          IF ( ln_slchldiainc ) THEN
+#if defined key_medusa
             zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
+#elif defined key_fabm
+            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1)
+#endif
             CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia )
          ELSE
@@ -954,4 +1320,30 @@
 #endif
 
+#if defined key_fabm
+         ! Nanophytoplankton chlorophyll
+         IF ( ln_slchlnaninc ) THEN
+            zbkg_chlnan(:,:) = tracer_bkg(:,:,1,jp_fabm_chl2)
+            CALL asm_bgc_unlog_2d( zbkg_chlnan, slchlnan_bkginc, zinc_chlnan )
+         ELSE
+            zinc_chlnan(:,:) = 0.0
+         ENDIF
+
+         ! Picophytoplankton chlorophyll
+         IF ( ln_slchlpicinc ) THEN
+            zbkg_chlpic(:,:) = tracer_bkg(:,:,1,jp_fabm_chl3)
+            CALL asm_bgc_unlog_2d( zbkg_chlpic, slchlpic_bkginc, zinc_chlpic )
+         ELSE
+            zinc_chlpic(:,:) = 0.0
+         ENDIF
+
+         ! Dinoflagellate chlorophyll
+         IF ( ln_slchldininc ) THEN
+            zbkg_chldin(:,:) = tracer_bkg(:,:,1,jp_fabm_chl4)
+            CALL asm_bgc_unlog_2d( zbkg_chldin, slchldin_bkginc, zinc_chldin )
+         ELSE
+            zinc_chldin(:,:) = 0.0
+         ENDIF
+#endif
+
          ! Total phytoplankton carbon
          IF ( ln_slphytotinc ) THEN
@@ -988,5 +1380,5 @@
          ! Select mixed layer
          IF ( ll_asmdin ) THEN
-#if defined key_top && ( defined key_hadocc || defined key_medusa )
+#if defined key_top && ( defined key_hadocc || defined key_medusa || defined key_fabm )
             CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
                &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
@@ -1014,7 +1406,5 @@
 #endif
             CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
-               !zmld(:,:) = hmld_tref(:,:)
-               CALL ctl_stop( ' hmld_tref mixed layer requested for phyto2d assimilation,', &
-                  &           ' but is not available in this version' )
+               zmld(:,:) = hmld_tref(:,:)
             CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
                zmld(:,:) = hmlpt(:,:)
@@ -1053,4 +1443,21 @@
             &                         cchl_p_bkg(:,:,1),                   &
             &                         tracer_bkg, phyto2d_balinc )
+#elif defined key_fabm
+         CALL asm_phyto2d_bal_ersem(  (ln_slchltotinc .OR. ln_schltotinc), &
+            &                         zinc_chltot,                         &
+            &                         ln_slchldiainc,                      &
+            &                         zinc_chldia,                         &
+            &                         ln_slchlnaninc,                      &
+            &                         zinc_chlnan,                         &
+            &                         ln_slchlpicinc,                      &
+            &                         zinc_chlpic,                         &
+            &                         ln_slchldininc,                      &
+            &                         zinc_chldin,                         &
+            &                         zincper,                             &
+            &                         rn_maxchlinc, ln_phytobal, zmld,     &
+            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
+            &                         phyt_avg_bkg, mld_max_bkg,           &
+            &                         totalk_bkg,                          &
+            &                         tracer_bkg, phyto2d_balinc )
 #else
          CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
@@ -1098,4 +1505,12 @@
                   &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
             END WHERE
+#elif defined key_fabm
+            WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
+               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
+            END WHERE
 #endif
 
@@ -1132,4 +1547,11 @@
                trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
             END WHERE
+#elif defined key_fabm
+            WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
+               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1)
+               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
+            END WHERE
 #endif
  
@@ -1163,14 +1585,31 @@
       REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
       !
-      INTEGER                          :: ji, jj, jk   ! Loop counters
-      INTEGER                          :: it           ! Index
-      REAL(wp)                         :: zincwgt      ! IAU weight for timestep
-      REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn
-      REAL(wp)                         :: zfrac_chd    ! Fraction of jpchd
-      REAL(wp)                         :: zrat_phn_chn ! jpphn:jpchn ratio
-      REAL(wp)                         :: zrat_phd_chd ! jpphd:jpchd ratio
-      REAL(wp)                         :: zrat_pds_chd ! jppds:jpchd ratio
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc      ! Chlorophyll increments
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl      ! Chlorophyll background
+      INTEGER                          :: ji, jj, jk    ! Loop counters
+      INTEGER                          :: it            ! Index
+      REAL(wp)                         :: zincwgt       ! IAU weight for timestep
+      REAL(wp)                         :: zfrac_chn     ! Fraction of jpchn
+      REAL(wp)                         :: zfrac_chd     ! Fraction of jpchd
+      REAL(wp)                         :: zfrac_chl1    ! Fraction of jp_fabm_chl1
+      REAL(wp)                         :: zfrac_chl2    ! Fraction of jp_fabm_chl2
+      REAL(wp)                         :: zfrac_chl3    ! Fraction of jp_fabm_chl3
+      REAL(wp)                         :: zfrac_chl4    ! Fraction of jp_fabm_chl4
+      REAL(wp)                         :: zrat_phn_chn  ! jpphn:jpchn ratio
+      REAL(wp)                         :: zrat_phd_chd  ! jpphd:jpchd ratio
+      REAL(wp)                         :: zrat_pds_chd  ! jppds:jpchd ratio
+      REAL(wp)                         :: zrat_p1c_chl1 ! jp_fabm_p1c:jp_fabm_chl1 ratio
+      REAL(wp)                         :: zrat_p1n_chl1 ! jp_fabm_p1n:jp_fabm_chl1 ratio
+      REAL(wp)                         :: zrat_p1p_chl1 ! jp_fabm_p1p:jp_fabm_chl1 ratio
+      REAL(wp)                         :: zrat_p1s_chl1 ! jp_fabm_p1s:jp_fabm_chl1 ratio
+      REAL(wp)                         :: zrat_p2c_chl2 ! jp_fabm_p2c:jp_fabm_chl2 ratio
+      REAL(wp)                         :: zrat_p2n_chl2 ! jp_fabm_p2n:jp_fabm_chl2 ratio
+      REAL(wp)                         :: zrat_p2p_chl2 ! jp_fabm_p2p:jp_fabm_chl2 ratio
+      REAL(wp)                         :: zrat_p3c_chl3 ! jp_fabm_p3c:jp_fabm_chl3 ratio
+      REAL(wp)                         :: zrat_p3n_chl3 ! jp_fabm_p3n:jp_fabm_chl3 ratio
+      REAL(wp)                         :: zrat_p3p_chl3 ! jp_fabm_p3p:jp_fabm_chl3 ratio
+      REAL(wp)                         :: zrat_p4c_chl4 ! jp_fabm_p4c:jp_fabm_chl4 ratio
+      REAL(wp)                         :: zrat_p4n_chl4 ! jp_fabm_p4n:jp_fabm_chl4 ratio
+      REAL(wp)                         :: zrat_p4p_chl4 ! jp_fabm_p4p:jp_fabm_chl4 ratio
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc       ! Chlorophyll increments
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl       ! Chlorophyll background
       !!------------------------------------------------------------------------
 
@@ -1189,4 +1628,9 @@
 #elif defined key_hadocc
             bkg_chl(:,:,:) = chl_bkg(:,:,:)
+#elif defined key_fabm
+            bkg_chl(:,:,:) = tracer_bkg(:,:,:,jp_fabm_chl1) + &
+               &             tracer_bkg(:,:,:,jp_fabm_chl2) + &
+               &             tracer_bkg(:,:,:,jp_fabm_chl3) + &
+               &             tracer_bkg(:,:,:,jp_fabm_chl4)
 #endif
             DO jk = 1, jpk
@@ -1240,4 +1684,51 @@
 #elif defined key_hadocc
          phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
+#elif defined key_fabm
+         ! Loop over each grid point partioning the increments based on existing ratios
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_chl1) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_chl2) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_chl3) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_chl4) > 0.0 ) ) THEN
+                     zfrac_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_chl1) / bkg_chl(ji,jj,jk)
+                     zfrac_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_chl2) / bkg_chl(ji,jj,jk)
+                     zfrac_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_chl3) / bkg_chl(ji,jj,jk)
+                     zfrac_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_chl4) / bkg_chl(ji,jj,jk)
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) = chl_inc(ji,jj,jk) * zfrac_chl1
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) = chl_inc(ji,jj,jk) * zfrac_chl2
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) = chl_inc(ji,jj,jk) * zfrac_chl3
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) = chl_inc(ji,jj,jk) * zfrac_chl4
+                     zrat_p1c_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1c) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
+                     zrat_p1n_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1n) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
+                     zrat_p1p_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1p) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
+                     zrat_p1s_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1s) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
+                     zrat_p2c_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2c) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
+                     zrat_p2n_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2n) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
+                     zrat_p2p_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2p) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
+                     zrat_p3c_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3c) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
+                     zrat_p3n_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3n) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
+                     zrat_p3p_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3p) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
+                     zrat_p4c_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4c) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
+                     zrat_p4n_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4n) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
+                     zrat_p4p_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4p) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p1c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1c_chl1
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p1n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1n_chl1
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p1p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1p_chl1
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p1s) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1s_chl1
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p2c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2c_chl2
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p2n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2n_chl2
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p2p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2p_chl2
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p3c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3c_chl3
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p3n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3n_chl3
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p3p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3p_chl3
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p4c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4c_chl4
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p4n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4n_chl4
+                     phyto3d_balinc(ji,jj,jk,jp_fabm_p4p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4p_chl4
+                  ENDIF
+               END DO
+            END DO
+         END DO
 #else
          CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
@@ -1285,4 +1776,12 @@
                   &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
             END WHERE
+#elif defined key_fabm
+            WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
+               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
+            END WHERE
 #endif
 
@@ -1318,4 +1817,11 @@
                   &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1)
                trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+            END WHERE
+#elif defined key_fabm
+            WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
+               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1)
+               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
             END WHERE
 #endif
@@ -1447,4 +1953,18 @@
             &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
 
+#elif defined key_fabm
+         ! Account for phytoplankton balancing if required
+         IF ( ln_phytobal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c) + phyto2d_balinc(:,:,1,jp_fabm_o3c)
+            alk_bkg_temp(:,:) = totalk_bkg(:,:,1)             + phyto2d_balinc(:,:,1,jp_fabm_o3ba)
+         ELSE
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c)
+            alk_bkg_temp(:,:) = totalk_bkg(:,:,1)
+         ENDIF
+
+         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
+            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
+            &               pco2_balinc(:,:,1,jp_fabm_o3c), pco2_balinc(:,:,1,jp_fabm_o3ba) )
+
 #else
          CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
@@ -1454,5 +1974,5 @@
          ! Select mixed layer
          IF ( ll_asmdin ) THEN
-#if defined key_hadocc || defined key_medusa
+#if defined key_hadocc || defined key_medusa || defined key_fabm
             CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
                &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
@@ -1480,7 +2000,5 @@
 #endif
             CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
-               !zmld(:,:) = hmld_tref(:,:)
-               CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
-                  &           ' but is not available in this version' )
+               zmld(:,:) = hmld_tref(:,:)
             CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
                zmld(:,:) = hmlpt(:,:)
@@ -1555,4 +2073,12 @@
                   &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
             END WHERE
+#elif defined key_fabm
+            WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
+               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
+            END WHERE
 #endif
 
@@ -1589,4 +2115,11 @@
                   &                         pco2_balinc(:,:,:,jp_had0:jp_had1)
                trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+            END WHERE
+#elif defined key_fabm
+            WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
+               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
+                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1)
+               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
             END WHERE
 #endif
@@ -1854,9 +2387,11 @@
 
          IF ( ln_pno3inc ) THEN
-#if defined key_hadocc || defined key_medusa
+#if defined key_hadocc || defined key_medusa || defined key_fabm
 #if defined key_hadocc
             it = jp_had_nut
 #elif defined key_medusa
             it = jpdin
+#elif defined key_fabm
+            it = jp_fabm_n3n
 #endif
             IF ( ln_phytobal ) THEN
@@ -1878,6 +2413,10 @@
 
          IF ( ln_psi4inc ) THEN
+#if defined key_medusa || defined key_fabm
 #if defined key_medusa
             it = jpsil
+#elif defined key_fabm
+            it = jp_fabm_n5s
+#endif
             IF ( ln_phytobal ) THEN
                psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
@@ -1898,9 +2437,11 @@
 
          IF ( ln_pdicinc ) THEN
-#if defined key_hadocc || defined key_medusa
+#if defined key_hadocc || defined key_medusa || defined key_fabm
 #if defined key_hadocc
             it = jp_had_dic
 #elif defined key_medusa
             it = jpdic
+#elif defined key_fabm
+            it = jp_fabm_o3c
 #endif
             IF ( ln_phytobal ) THEN
@@ -1922,9 +2463,11 @@
 
          IF ( ln_palkinc ) THEN
-#if defined key_hadocc || defined key_medusa
+#if defined key_hadocc || defined key_medusa || defined key_fabm
 #if defined key_hadocc
             it = jp_had_alk
 #elif defined key_medusa
             it = jpalk
+#elif defined key_fabm
+            it = jp_fabm_o3ba
 #endif
             IF ( ln_phytobal ) THEN
@@ -1946,6 +2489,10 @@
 
          IF ( ln_po2inc ) THEN
+#if defined key_medusa || defined key_fabm
 #if defined key_medusa
             it = jpoxy
+#elif defined key_fabm
+            it = jp_fabm_o2o
+#endif
             IF ( ln_phytobal ) THEN
                po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
@@ -2004,4 +2551,10 @@
                   trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
                END WHERE
+#elif defined key_fabm
+               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jp_fabm_n3n) = trb(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt
+               END WHERE
 #else
                CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
@@ -2015,4 +2568,10 @@
                   trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
                   trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
+               END WHERE
+#elif defined key_fabm
+               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jp_fabm_n5s) = trb(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt
                END WHERE
 #else
@@ -2034,4 +2593,10 @@
                   trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
                END WHERE
+#elif defined key_fabm
+               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jp_fabm_o3c) = trb(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt
+               END WHERE
 #else
                CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
@@ -2052,4 +2617,10 @@
                   trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
                END WHERE
+#elif defined key_fabm
+               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jp_fabm_o3ba) = trb(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt
+               END WHERE
 #else
                CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
@@ -2063,4 +2634,10 @@
                   trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
                   trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
+               END WHERE
+#elif defined key_fabm
+               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jp_fabm_o2o) = trb(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt
                END WHERE
 #else
@@ -2091,11 +2668,6 @@
             ! Initialize the now fields with the background + increment
             ! Background currently is what the model is initialised with
-#if defined key_hadocc
-            CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', &
+            CALL ctl_warn( ' Doing direct initialisation with 3D BGC assimilation', &
                &           ' Background state is taken from model rather than background file' )
-#elif defined key_medusa
-            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', &
-               &           ' Background state is taken from model rather than background file' )
-#endif
 
             IF ( ln_pno3inc ) THEN
@@ -2112,4 +2684,10 @@
                   trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
                END WHERE
+#elif defined key_fabm
+               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:)
+                  trb(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n)
+               END WHERE
 #else
                CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
@@ -2123,4 +2701,10 @@
                   trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
                   trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
+               END WHERE
+#elif defined key_fabm
+               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:)
+                  trb(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s)
                END WHERE
 #else
@@ -2142,4 +2726,10 @@
                   trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
                END WHERE
+#elif defined key_fabm
+               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:)
+                  trb(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c)
+               END WHERE
 #else
                CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
@@ -2160,4 +2750,10 @@
                   trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
                END WHERE
+#elif defined key_fabm
+               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:)
+                  trb(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba)
+               END WHERE
 #else
                CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
@@ -2171,4 +2767,10 @@
                   trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
                   trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
+               END WHERE
+#elif defined key_fabm
+               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:)
+                  trb(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o)
                END WHERE
 #else
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 10574)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 10622)
@@ -52,5 +52,7 @@
    USE asminc, ONLY: ln_avgbkg
 #if defined key_top
-   USE asmbgc, ONLY: asm_bgc_bkg_wri
+   USE asmbgc, ONLY: asm_bgc_bkg_alloc, &
+      &              asm_bgc_bkg_tavg,  &
+      &              asm_bgc_bkg_wri
 #endif
    IMPLICIT NONE
@@ -140,5 +142,12 @@
          
          numtimes_tavg = REAL ( nitavgbkg_r -  nn_it000 + 1 )
-      ENDIF   
+      ENDIF
+      
+#if defined key_top
+      ! Allocate BGC average arrays whatever, to save code repetition later
+      IF ( kt == ( nn_it000 - 1) ) THEN
+         CALL asm_bgc_bkg_alloc
+      ENDIF
+#endif
 
       ! If creating an averaged assim bkg, sum the contribution every timestep
@@ -157,4 +166,7 @@
 #if defined key_zdftke
          en_tavg(:,:,:)       = en_tavg(:,:,:) + en(:,:,:) / numtimes_tavg
+#endif
+#if defined key_top
+         CALL asm_bgc_bkg_tavg( kt, numtimes_tavg )
 #endif
       ENDIF
@@ -226,5 +238,5 @@
             
 #if defined key_top
-            CALL asm_bgc_bkg_wri( kt, inum )
+            CALL asm_bgc_bkg_wri( kt, inum, ln_avgbkg )
 #endif
             CALL iom_close( inum )
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 10574)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 10622)
@@ -167,4 +167,5 @@
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
          &                 ln_phytobal, ln_slchltotinc, ln_slchldiainc,    &
+         &                 ln_slchlnaninc, ln_slchlpicinc, ln_slchldininc, &
          &                 ln_slchlnoninc, ln_schltotinc, ln_slphytotinc,  &
          &                 ln_slphydiainc, ln_slphynoninc, ln_spco2inc,    &
@@ -228,4 +229,7 @@
          WRITE(numout,*) '      Logical switch for applying slchldia increments     ln_slchldiainc = ', ln_slchldiainc
          WRITE(numout,*) '      Logical switch for applying slchlnon increments     ln_slchlnoninc = ', ln_slchlnoninc
+         WRITE(numout,*) '      Logical switch for applying slchlnan increments     ln_slchlnaninc = ', ln_slchlnaninc
+         WRITE(numout,*) '      Logical switch for applying slchlpic increments     ln_slchlpicinc = ', ln_slchlpicinc
+         WRITE(numout,*) '      Logical switch for applying slchldin increments     ln_slchldininc = ', ln_slchldininc
          WRITE(numout,*) '      Logical switch for applying schltot increments       ln_schltotinc = ', ln_schltotinc
          WRITE(numout,*) '      Logical switch for applying slphytot increments     ln_slphytotinc = ', ln_slphytotinc
@@ -295,4 +299,5 @@
       ENDIF
       IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
          & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
          & ln_slphynoninc .OR. ln_spco2inc    .OR. ln_sfco2inc    .OR. &
@@ -1353,4 +1358,5 @@
       ! Ocean colour variables first
       IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
          & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
          & ln_slphynoninc ) THEN
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_ersem.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_ersem.F90	(revision 10622)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_ersem.F90	(revision 10622)
@@ -0,0 +1,747 @@
+MODULE asmphyto2dbal_ersem
+   !!======================================================================
+   !!                       ***  MODULE asmphyto2dbal_ersem  ***
+   !! Calculate increments to ERSEM based on surface phyto2d increments
+   !!
+   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
+   !! For licensing reasons this is kept in its own internal Met Office
+   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
+   !! repository, and must be merged in when building.
+   !!
+   !!======================================================================
+   !! History :  3.6  ! 2019-01 (D. Ford)  Adapted from asmphyto2dbal_medusa
+   !!----------------------------------------------------------------------
+#if defined key_asminc && defined key_fabm
+   !!----------------------------------------------------------------------
+   !! 'key_asminc'          : assimilation increment interface
+   !! 'key_fabm'            : FABM-ERSEM model
+   !!----------------------------------------------------------------------
+   !! asm_phyto2d_bal_ersem : routine to calculate increments to ERSEM
+   !!----------------------------------------------------------------------
+   USE par_kind,      ONLY: wp             ! kind parameters
+   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
+   USE dom_oce,       ONLY: gdepw_n        ! domain information
+   USE iom                                 ! i/o
+   USE par_fabm                            ! FABM-ERSEM parameters
+   USE par_trc,       ONLY: jptra          ! Tracer parameters
+   USE bioanalysis                         ! Nitrogen balancing
+
+   IMPLICIT NONE
+   PRIVATE                   
+
+   PUBLIC asm_phyto2d_bal_ersem
+
+   ! Default values for biological assimilation parameters
+   ! Should match Hemmings et al. (2008)
+   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
+   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
+   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
+   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
+   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
+   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
+   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
+   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
+   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
+   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
+   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
+   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
+   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
+   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
+   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
+   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
+   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
+   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
+   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
+   !
+   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
+   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
+   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
+   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
+   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
+   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
+   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
+   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
+
+CONTAINS
+
+   SUBROUTINE asm_phyto2d_bal_ersem(  ld_chltot,                      &
+      &                               pinc_chltot,                    &
+      &                               ld_chldia,                      &
+      &                               pinc_chldia,                    &
+      &                               ld_chlnan,                      &
+      &                               pinc_chlnan,                    &
+      &                               ld_chlpic,                      &
+      &                               pinc_chlpic,                    &
+      &                               ld_chldin,                      &
+      &                               pinc_chldin,                    &
+      &                               pincper,                        &
+      &                               p_maxchlinc, ld_phytobal, pmld, &
+      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                               phyt_avg_bkg, mld_max_bkg,      &
+      &                               totalk_bkg,                     &
+      &                               tracer_bkg, phyto2d_balinc )
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_phyto2d_bal_ersem  ***
+      !!
+      !! ** Purpose :   calculate increments to ERSEM from 2d phytoplankton increments
+      !!
+      !! ** Method  :   EITHER (ld_phytobal == .TRUE.):
+      !!                average up ERSEM to look like HadOCC
+      !!                call nitrogen balancing scheme
+      !!                separate back out to MEDUSA
+      !!                OR (ld_phytobal == .FALSE.):
+      !!                calculate increments to maintain background stoichiometry
+      !!
+      !! ** Action  :   populate phyto2d_balinc
+      !!
+      !! References :   Hemmings et al., 2008, J. Mar. Res.
+      !!                Ford et al., 2012, Ocean Sci.
+      !!                Skakala et al., 2018, JGR
+      !!---------------------------------------------------------------------------
+      !!
+      LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments
+      LOGICAL,  INTENT(in   )                               :: ld_chldia      ! Assim chldia y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldia    ! chldia increments
+      LOGICAL,  INTENT(in   )                               :: ld_chlnan      ! Assim chlnan y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlnan    ! chlnan increments
+      LOGICAL,  INTENT(in   )                               :: ld_chlpic      ! Assim chlpic y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlpic    ! chlpic increments
+      LOGICAL,  INTENT(in   )                               :: ld_chldin      ! Assim chldin y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldin    ! chldin increments
+      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
+      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
+      LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: totalk_bkg     ! Total alkalinity
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
+      !!
+      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
+      INTEGER                                               :: jkmax          ! Loop index
+      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
+      REAL(wp)                                              :: zmassc         ! Carbon molar mass
+      REAL(wp)                                              :: zmassn         ! Nitrogen molar mass
+      REAL(wp)                                              :: z4qnc          ! Z4/qnc (mesozoo N:C)
+      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
+      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
+      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
+      REAL(wp)                                              :: zfrac          ! Fractions
+      REAL(wp)                                              :: zfrac_chl1     ! 
+      REAL(wp)                                              :: zfrac_chl2     ! 
+      REAL(wp)                                              :: zfrac_chl3     ! 
+      REAL(wp)                                              :: zfrac_chl4     ! 
+      REAL(wp)                                              :: zfrac_p1n      ! 
+      REAL(wp)                                              :: zfrac_p2n      ! 
+      REAL(wp)                                              :: zfrac_p3n      ! 
+      REAL(wp)                                              :: zfrac_p4n      ! 
+      REAL(wp)                                              :: zfrac_z4n      ! 
+      REAL(wp)                                              :: zfrac_z5n      ! 
+      REAL(wp)                                              :: zfrac_z6n      ! 
+      REAL(wp)                                              :: zfrac_n3n      ! 
+      REAL(wp)                                              :: zfrac_n4n      ! 
+      REAL(wp)                                              :: zfrac_r4n      ! 
+      REAL(wp)                                              :: zfrac_r6n      ! 
+      REAL(wp)                                              :: zfrac_r8n      ! 
+      REAL(wp)                                              :: zrat_chl1_p1n  ! Ratios
+      REAL(wp)                                              :: zrat_p1c_p1n   ! 
+      REAL(wp)                                              :: zrat_p1p_p1n   ! 
+      REAL(wp)                                              :: zrat_p1s_p1n   ! 
+      REAL(wp)                                              :: zrat_chl2_p2n  ! 
+      REAL(wp)                                              :: zrat_p2c_p2n   ! 
+      REAL(wp)                                              :: zrat_p2p_p2n   ! 
+      REAL(wp)                                              :: zrat_chl3_p3n  ! 
+      REAL(wp)                                              :: zrat_p3c_p3n   ! 
+      REAL(wp)                                              :: zrat_p3p_p3n   ! 
+      REAL(wp)                                              :: zrat_chl4_p4n  ! 
+      REAL(wp)                                              :: zrat_p4c_p4n   ! 
+      REAL(wp)                                              :: zrat_p4p_p4n   ! 
+      REAL(wp)                                              :: zrat_z4c_z4n   ! 
+      REAL(wp)                                              :: zrat_z5c_z5n   ! 
+      REAL(wp)                                              :: zrat_z5p_z5n   ! 
+      REAL(wp)                                              :: zrat_z6c_z6n   ! 
+      REAL(wp)                                              :: zrat_z6p_z6n   ! 
+      REAL(wp)                                              :: zrat_r4c_r4n   ! 
+      REAL(wp)                                              :: zrat_r4p_r4n   ! 
+      REAL(wp)                                              :: zrat_r6c_r6n   ! 
+      REAL(wp)                                              :: zrat_r6p_r6n   ! 
+      REAL(wp)                                              :: zrat_r6s_r6n   ! 
+      REAL(wp)                                              :: zrat_r8c_r8n   ! 
+      REAL(wp)                                              :: zrat_r8p_r8n   ! 
+      REAL(wp)                                              :: zrat_r8s_r8n   ! 
+      REAL(wp)                                              :: zrat_p1c_chl1  ! 
+      REAL(wp)                                              :: zrat_p1n_chl1  ! 
+      REAL(wp)                                              :: zrat_p1p_chl1  ! 
+      REAL(wp)                                              :: zrat_p1s_chl1  ! 
+      REAL(wp)                                              :: zrat_p2c_chl2  ! 
+      REAL(wp)                                              :: zrat_p2n_chl2  ! 
+      REAL(wp)                                              :: zrat_p2p_chl2  ! 
+      REAL(wp)                                              :: zrat_p3c_chl3  ! 
+      REAL(wp)                                              :: zrat_p3n_chl3  ! 
+      REAL(wp)                                              :: zrat_p3p_chl3  ! 
+      REAL(wp)                                              :: zrat_p4c_chl4  ! 
+      REAL(wp)                                              :: zrat_p4n_chl4  ! 
+      REAL(wp)                                              :: zrat_p4p_chl4  ! 
+      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
+      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
+      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
+      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_temp
+      !!---------------------------------------------------------------------------
+      
+      ! Set parameters
+      zmassc = 12.01
+      zmassn = 14.01
+      z4qnc  = 0.0126
+      !z4qnc  = model%state_variables(jp_fabm_z4c)%parameters%qnc%value
+      !z4qnc  = get_property_by_name(model%state_variables(jp_fabm_z4c)%parameters, 'qnc')
+      IF (lwp) WRITE(numout,*) 'z4qnc = ', z4qnc
+
+      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      IF ( p_maxchlinc > 0.0 ) THEN
+         IF ( ld_chltot ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
+               END DO
+            END DO
+         ELSE
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
+                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
+                        &                      pinc_chlpic(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
+                        &                      pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlpic .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + &
+                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + &
+                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlnan ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlpic ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlpic(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chlnan .AND. ld_chlpic ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chlpic(ji,jj)
+                  ELSE IF ( ld_chlnan .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chlpic .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chldia ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj)
+                  ELSE IF ( ld_chlnan ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj)
+                  ELSE IF ( ld_chlpic ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj)
+                  ELSE IF ( ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldin(ji,jj)
+                  ENDIF
+                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_temp(ji,jj), p_maxchlinc ) )
+                  IF ( pinc_chltot(ji,jj) .NE. pinc_chltot_temp(ji,jj) ) THEN
+                     zfrac = pinc_chltot(ji,jj) / pinc_chltot_temp(ji,jj)
+                     IF ( ld_chldia ) THEN
+                        pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac
+                     ENDIF
+                     IF ( ld_chlnan ) THEN
+                        pinc_chlnan(ji,jj) = pinc_chlnan(ji,jj) * zfrac
+                     ENDIF
+                     IF ( ld_chlpic ) THEN
+                        pinc_chlpic(ji,jj) = pinc_chlpic(ji,jj) * zfrac
+                     ENDIF
+                     IF ( ld_chldin ) THEN
+                        pinc_chldin(ji,jj) = pinc_chldin(ji,jj) * zfrac
+                     ENDIF
+                  ENDIF
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+      
+      ! Initialise balancing increments
+      phyto2d_balinc(:,:,:,:) = 0.0
+      
+      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
+
+         ! Set up model parameters to be passed into Hemmings balancing routine.
+         ! For now these are hardwired to the standard HadOCC parameter values
+         ! as this is what the scheme was developed for.
+         ! Obviously, HadOCC and ERSEM are rather different models, so this
+         ! isn't ideal, but there's not always direct analogues between the two
+         ! parameter sets, so it's the easiest way to get something running.
+         ! In the longer term, some serious MarMOT-based development is required.
+         modparm(1)  = 0.1                               ! grow_sat
+         modparm(2)  = 2.0                               ! psmax
+         modparm(3)  = 0.845                             ! par
+         modparm(4)  = 0.02                              ! alpha
+         modparm(5)  = 0.05                              ! resp_rate
+         modparm(6)  = 0.05                              ! pmort_rate
+         modparm(7)  = 0.01                              ! phyto_min
+         modparm(8)  = 0.05                              ! z_mort_1
+         modparm(9)  = 1.0                               ! z_mort_2
+         modparm(10) = 6.625                             ! c2n_p
+         modparm(11) = 5.625                             ! c2n_z
+         modparm(12) = 7.5                               ! c2n_d
+         modparm(13) = 0.01                              ! graze_threshold
+         modparm(14) = 2.0                               ! holling_coef
+         modparm(15) = 0.5                               ! graze_sat
+         modparm(16) = 2.0                               ! graze_max
+
+         ! Set up assimilation parameters to be passed into balancing routine
+         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
+         assimparm(2)  = balnutext
+         assimparm(3)  = balnutmin
+         assimparm(4)  = r
+         assimparm(5)  = beta_g
+         assimparm(6)  = beta_l
+         assimparm(7)  = beta_m
+         assimparm(8)  = a_g
+         assimparm(9)  = a_l
+         assimparm(10) = a_m
+         assimparm(11) = zfracb0
+         assimparm(12) = zfracb1
+         assimparm(13) = qrfmax
+         assimparm(14) = qafmax
+         assimparm(15) = zrfmax
+         assimparm(16) = zafmax
+         assimparm(17) = prfmax
+         assimparm(18) = incphymin
+         assimparm(19) = integnstep
+         assimparm(20) = pthreshold
+
+         ! Set up external tracer indices array bstate
+         i_tracer(1) = 1   ! nutrient
+         i_tracer(2) = 2   ! phytoplankton
+         i_tracer(3) = 3   ! zooplankton
+         i_tracer(4) = 4   ! detritus
+         i_tracer(5) = 5   ! DIC
+         i_tracer(6) = 6   ! Alkalinity
+
+         ! Set background state
+         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_fabm_n3n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_n4n)
+         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_fabm_p1n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_p2n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_p3n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_p4n)
+         bstate(:,:,:,i_tracer(3)) = (tracer_bkg(:,:,:,jp_fabm_z4c) * z4qnc) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_z5n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_z6n)
+         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_fabm_r4n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_r6n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_r8n)
+         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_fabm_o3c)
+         bstate(:,:,:,i_tracer(6)) = totalk_bkg(:,:,:)
+
+         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
+         ! and nitrogen to biomass equivalent for PZD (hardwire as per HadOCC)
+         cchl_p(:,:) = 0.0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                  &   tracer_bkg(ji,jj,1,jp_fabm_chl3) + tracer_bkg(ji,jj,1,jp_fabm_chl4) ) .GT. 0.0 ) THEN
+                  cchl_p(ji,jj) = zmassc * ( tracer_bkg(ji,jj,1,jp_fabm_p1c) +      &
+                     &                       tracer_bkg(ji,jj,1,jp_fabm_p2c) +      &
+                     &                       tracer_bkg(ji,jj,1,jp_fabm_p3c) +      &
+                     &                       tracer_bkg(ji,jj,1,jp_fabm_p4c)   ) /  &
+                     &            ( tracer_bkg(ji,jj,1,jp_fabm_chl1) +             &
+                     &              tracer_bkg(ji,jj,1,jp_fabm_chl2) +             &
+                     &              tracer_bkg(ji,jj,1,jp_fabm_chl3) +             &
+                     &              tracer_bkg(ji,jj,1,jp_fabm_chl4)   )
+               ENDIF
+            END DO
+         END DO
+         n2be_p = zmassn + ( zmassc * modparm(10) )
+         n2be_z = zmassn + ( zmassc * modparm(11) )
+         n2be_d = zmassn + ( zmassc * modparm(12) )
+
+         ! Call nitrogen balancing routine
+         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
+            &               n2be_p, n2be_z, n2be_d, assimparm,                      &
+            &               INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
+            &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p(:,:), &
+            &               nbal_active, phyt_avg_bkg(:,:),                         &
+            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      &
+            &               subsurf_active, deepneg_active,                         &
+            &               deeppos_active, nutprof_active,                         &
+            &               bstate, outincs,                                        &
+            &               diag_active, diag,                                      &
+            &               diag_fulldepth_active, diag_fulldepth )
+         
+         ! Loop over each grid point partioning the increments
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  ! Phytoplankton
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p2n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p3n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p4n) > 0.0 ) .AND. &
+                     & ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN
+                     IF ( ld_chltot ) THEN
+                        ! Phytoplankton nitrogen split up based on existing ratios
+                        zfrac_p1n = tracer_bkg(ji,jj,jk,jp_fabm_p1n) / &
+                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) )
+                        zfrac_p2n = tracer_bkg(ji,jj,jk,jp_fabm_p2n) / &
+                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) )
+                        zfrac_p3n = tracer_bkg(ji,jj,jk,jp_fabm_p3n) / &
+                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) )
+                        zfrac_p4n = tracer_bkg(ji,jj,jk,jp_fabm_p4n) / &
+                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) )
+                     ELSE
+                        ! Phytoplankton nitrogen split up based on assimilation increments
+                        zfrac_p1n = pinc_chldia(ji,jj) / pinc_chltot(ji,jj)
+                        zfrac_p2n = pinc_chlnan(ji,jj) / pinc_chltot(ji,jj)
+                        zfrac_p3n = pinc_chlpic(ji,jj) / pinc_chltot(ji,jj)
+                        zfrac_p4n = pinc_chldin(ji,jj) / pinc_chltot(ji,jj)
+                     ENDIF
+                     
+                     ! Other phytoplankton variables split up based on existing ratios with nitrogen
+                     zrat_chl1_p1n = tracer_bkg(ji,jj,jk,jp_fabm_chl1) / tracer_bkg(ji,jj,jk,jp_fabm_p1n)
+                     zrat_p1c_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1c)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n)
+                     zrat_p1p_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1p)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n)
+                     zrat_p1s_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1s)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n)
+                     zrat_chl2_p2n = tracer_bkg(ji,jj,jk,jp_fabm_chl2) / tracer_bkg(ji,jj,jk,jp_fabm_p2n)
+                     zrat_p2c_p2n  = tracer_bkg(ji,jj,jk,jp_fabm_p2c)  / tracer_bkg(ji,jj,jk,jp_fabm_p2n)
+                     zrat_p2p_p2n  = tracer_bkg(ji,jj,jk,jp_fabm_p2p)  / tracer_bkg(ji,jj,jk,jp_fabm_p2n)
+                     zrat_chl3_p3n = tracer_bkg(ji,jj,jk,jp_fabm_chl3) / tracer_bkg(ji,jj,jk,jp_fabm_p3n)
+                     zrat_p3c_p3n  = tracer_bkg(ji,jj,jk,jp_fabm_p3c)  / tracer_bkg(ji,jj,jk,jp_fabm_p3n)
+                     zrat_p3p_p3n  = tracer_bkg(ji,jj,jk,jp_fabm_p3p)  / tracer_bkg(ji,jj,jk,jp_fabm_p3n)
+                     zrat_chl4_p4n = tracer_bkg(ji,jj,jk,jp_fabm_chl4) / tracer_bkg(ji,jj,jk,jp_fabm_p4n)
+                     zrat_p4c_p4n  = tracer_bkg(ji,jj,jk,jp_fabm_p4c)  / tracer_bkg(ji,jj,jk,jp_fabm_p4n)
+                     zrat_p4p_p4n  = tracer_bkg(ji,jj,jk,jp_fabm_p4p)  / tracer_bkg(ji,jj,jk,jp_fabm_p4n)
+                     
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl1) = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_chl1_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1c_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1p_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1s)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1s_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl2) = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_chl2_p2n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_p2c_p2n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_p2p_p2n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl3) = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_chl3_p3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_p3c_p3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_p3p_p3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl4) = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_chl4_p4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_p4c_p4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_p4p_p4n
+                  ENDIF
+
+                  ! Zooplankton nitrogen split up based on existing ratios
+                  ! Update carbon and phosphorus according to existing ratios
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_z5n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_z6n) > 0.0 ) ) THEN
+                     zfrac_z4n = ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) / &
+                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) )
+                     zfrac_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5n) / &
+                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) )
+                     zfrac_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6n) / &
+                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) )
+                     zrat_z4c_z4n = 1.0 / z4qnc
+                     zrat_z5c_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5c) / tracer_bkg(ji,jj,jk,jp_fabm_z5n)
+                     zrat_z5p_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5p) / tracer_bkg(ji,jj,jk,jp_fabm_z5n)
+                     zrat_z6c_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6c) / tracer_bkg(ji,jj,jk,jp_fabm_z6n)
+                     zrat_z6p_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6p) / tracer_bkg(ji,jj,jk,jp_fabm_z6n)
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5c) = phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) * zrat_z5c_z5n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6c) = phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) * zrat_z6c_z6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5p) = phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) * zrat_z5p_z5n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6p) = phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) * zrat_z6p_z6n
+                  ENDIF
+
+                  ! Nitrogen nutrient split between nitrate and ammonium based on existing ratios
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n3n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_n4n) > 0.0 ) ) THEN
+                     zfrac_n3n = tracer_bkg(ji,jj,jk,jp_fabm_n3n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_n4n))
+                     zfrac_n4n = tracer_bkg(ji,jj,jk,jp_fabm_n4n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_n4n))
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n
+                  ENDIF
+
+                  ! Detritus nitrogen split up based on existing ratios
+                  ! Update carbon and phosphorus according to existing ratios
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_r4n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_r6n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_r8n) > 0.0 ) ) THEN
+                     zfrac_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n))
+                     zfrac_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n))
+                     zfrac_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n))
+                     zrat_r4c_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4c) / tracer_bkg(ji,jj,jk,jp_fabm_r4n)
+                     zrat_r4p_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4p) / tracer_bkg(ji,jj,jk,jp_fabm_r4n)
+                     zrat_r6c_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6c) / tracer_bkg(ji,jj,jk,jp_fabm_r6n)
+                     zrat_r6p_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6p) / tracer_bkg(ji,jj,jk,jp_fabm_r6n)
+                     zrat_r6s_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6s) / tracer_bkg(ji,jj,jk,jp_fabm_r6n)
+                     zrat_r8c_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8c) / tracer_bkg(ji,jj,jk,jp_fabm_r8n)
+                     zrat_r8p_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8p) / tracer_bkg(ji,jj,jk,jp_fabm_r8n)
+                     zrat_r8s_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8s) / tracer_bkg(ji,jj,jk,jp_fabm_r8n)
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) * zrat_r4c_r4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) * zrat_r4p_r4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6c_r6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6p_r6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6s) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6s_r6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8c_r8n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8p_r8n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8s) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8s_r8n
+                  ENDIF
+
+                  ! DIC straight from balancing scheme
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5))
+
+                  ! Alkalinity straight from balancing scheme
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_o3ba) = outincs(ji,jj,jk,i_tracer(6))
+
+                  ! Remove P/R silicon increments from silicate to conserve mass
+                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_p1s) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r6s) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r8s)
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n5s) - zfrac ) > 0.0 ) THEN
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_n5s) = zfrac * (-1.0)
+                  ENDIF
+
+                  ! Remove P/Z/R phosphorus increments from phosphate to conserve mass
+                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_p1p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p2p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p3p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p4p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_z5p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_z6p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r4p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r6p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r8p)
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n1p) - zfrac ) > 0.0 ) THEN
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_n1p) = zfrac * (-1.0)
+                  ENDIF
+                  
+               END DO
+            END DO
+         END DO
+      
+      ELSE   ! No nitrogen balancing - just update phytoplankton
+         
+         ! Split up total surface chlorophyll increments
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_chl1) > 0.0 ) .AND. &
+                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl2) > 0.0 ) .AND. &
+                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl3) > 0.0 ) .AND. &
+                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl4) > 0.0 ) ) THEN
+                  IF ( ld_chltot ) THEN
+                     ! Chlorophyll split up based on existing ratios
+                     zfrac_chl1 = tracer_bkg(ji,jj,1,jp_fabm_chl1) /  &
+                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) )
+                     zfrac_chl2 = tracer_bkg(ji,jj,1,jp_fabm_chl2) /  &
+                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) )
+                     zfrac_chl3 = tracer_bkg(ji,jj,1,jp_fabm_chl3) /  &
+                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) )
+                     zfrac_chl4 = tracer_bkg(ji,jj,1,jp_fabm_chl4) /  &
+                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) )
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl1) = pinc_chltot(ji,jj) * zfrac_chl1
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl2) = pinc_chltot(ji,jj) * zfrac_chl2
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl3) = pinc_chltot(ji,jj) * zfrac_chl3
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl4) = pinc_chltot(ji,jj) * zfrac_chl4
+                  ENDIF
+                  IF( ld_chldia ) THEN
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl1) = pinc_chldia(ji,jj)
+                  ENDIF
+                  IF( ld_chlnan ) THEN
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl2) = pinc_chlnan(ji,jj)
+                  ENDIF
+                  IF( ld_chlpic ) THEN
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl3) = pinc_chlpic(ji,jj)
+                  ENDIF
+                  IF( ld_chldin ) THEN
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl4) = pinc_chldin(ji,jj)
+                  ENDIF
+                  
+                  ! Maintain stoichiometric ratios of carbon, nitrogen, phosphorus and silicon
+                  IF ( ld_chltot .OR. ld_chldia ) THEN
+                     zrat_p1c_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1c) / tracer_bkg(ji,jj,1,jp_fabm_chl1)
+                     zrat_p1n_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1n) / tracer_bkg(ji,jj,1,jp_fabm_chl1)
+                     zrat_p1p_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1p) / tracer_bkg(ji,jj,1,jp_fabm_chl1)
+                     zrat_p1s_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1s) / tracer_bkg(ji,jj,1,jp_fabm_chl1)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p1c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1c_chl1
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p1n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1n_chl1
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p1p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1p_chl1
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p1s) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1s_chl1
+                  ENDIF
+                  IF ( ld_chltot .OR. ld_chlnan ) THEN
+                     zrat_p2c_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2c) / tracer_bkg(ji,jj,1,jp_fabm_chl2)
+                     zrat_p2n_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2n) / tracer_bkg(ji,jj,1,jp_fabm_chl2)
+                     zrat_p2p_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2p) / tracer_bkg(ji,jj,1,jp_fabm_chl2)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p2c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2c_chl2
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p2n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2n_chl2
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p2p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2p_chl2
+                  ENDIF
+                  IF ( ld_chltot .OR. ld_chlpic ) THEN
+                     zrat_p3c_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3c) / tracer_bkg(ji,jj,1,jp_fabm_chl3)
+                     zrat_p3n_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3n) / tracer_bkg(ji,jj,1,jp_fabm_chl3)
+                     zrat_p3p_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3p) / tracer_bkg(ji,jj,1,jp_fabm_chl3)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p3c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3c_chl3
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p3n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3n_chl3
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p3p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3p_chl3
+                  ENDIF
+                  IF ( ld_chltot .OR. ld_chldin ) THEN
+                     zrat_p4c_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4c) / tracer_bkg(ji,jj,1,jp_fabm_chl4)
+                     zrat_p4n_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4n) / tracer_bkg(ji,jj,1,jp_fabm_chl4)
+                     zrat_p4p_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4p) / tracer_bkg(ji,jj,1,jp_fabm_chl4)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p4c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4c_chl4
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p4n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4n_chl4
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p4p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4p_chl4
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+         
+         ! Propagate through mixed layer
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               jkmax = jpk-1
+               DO jk = jpk-1, 1, -1
+                  IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                     & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                     pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                     jkmax = jk
+                  ENDIF
+               END DO
+               !
+               DO jk = 2, jkmax
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl1) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1c)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1n)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1p)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1s)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1s)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl2) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2c)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2n)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2p)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl3) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3c)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3n)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3p)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl4) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4c)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4n)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4p)
+               END DO
+               !
+            END DO
+         END DO
+
+      ENDIF
+
+   END SUBROUTINE asm_phyto2d_bal_ersem
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option : Empty routine
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE asm_phyto2d_bal_ersem(  ld_chltot,                      &
+      &                               pinc_chltot,                    &
+      &                               ld_chldia,                      &
+      &                               pinc_chldia,                    &
+      &                               ld_chlnan,                      &
+      &                               pinc_chlnan,                    &
+      &                               ld_chlpic,                      &
+      &                               pinc_chlpic,                    &
+      &                               ld_chldin,                      &
+      &                               pinc_chldin,                    &
+      &                               pincper,                        &
+      &                               p_maxchlinc, ld_phytobal, pmld, &
+      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                               phyt_avg_bkg, mld_max_bkg,      &
+      &                               totalk_bkg,                     &
+      &                               tracer_bkg, phyto2d_balinc )
+      LOGICAL :: ld_chltot
+      REAL    :: pinc_chltot(:,:)
+      LOGICAL :: ld_chldia
+      REAL    :: pinc_chldia(:,:)
+      LOGICAL :: ld_chlnan
+      REAL    :: pinc_chlnan(:,:)
+      LOGICAL :: ld_chlpic
+      REAL    :: pinc_chlpic(:,:)
+      LOGICAL :: ld_chldin
+      REAL    :: pinc_chldin(:,:)
+      REAL    :: pincper
+      REAL    :: p_maxchlinc
+      LOGICAL :: ld_phytobal
+      REAL    :: pmld(:,:)
+      REAL    :: pgrow_avg_bkg(:,:)
+      REAL    :: ploss_avg_bkg(:,:)
+      REAL    :: phyt_avg_bkg(:,:)
+      REAL    :: mld_max_bkg(:,:)
+      REAL    :: totalk_bkg(:,:,:)
+      REAL    :: tracer_bkg(:,:,:,:)
+      REAL    :: phyto2d_balinc(:,:,:,:)
+      WRITE(*,*) 'asm_phyto2d_bal_ersem: You should not have seen this print! error?'
+   END SUBROUTINE asm_phyto2d_bal_ersem
+#endif
+
+   !!======================================================================
+END MODULE asmphyto2dbal_ersem
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 10574)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90	(revision 10622)
@@ -12,21 +12,28 @@
 
    ! Variables needed for OBS/ASM
-   INTEGER, PUBLIC :: jp_fabm_chl1, jp_fabm_chl2, &
-                      jp_fabm_chl3, jp_fabm_chl4, &
-                      jp_fabm_p1c,  jp_fabm_p1n,  &
-                      jp_fabm_p1p,  jp_fabm_p1s,  &
-                      jp_fabm_p2c,  jp_fabm_p2n,  &
-                      jp_fabm_p2p,  jp_fabm_p3c,  &
-                      jp_fabm_p3n,  jp_fabm_p3p,  &
-                      jp_fabm_p4c,  jp_fabm_p4n,  &
-                      jp_fabm_p4p,  jp_fabm_z4c,  &
-                      jp_fabm_z5c,  jp_fabm_z5n,  &
-                      jp_fabm_z5p,  jp_fabm_z6c,  &
-                      jp_fabm_z6n,  jp_fabm_z6p,  &
-                      jp_fabm_n1p,  jp_fabm_n3n,  &
-                      jp_fabm_n4n,  jp_fabm_n5s,  &
-                      jp_fabm_o2o,  jp_fabm_o3c,  &
-                      jp_fabm_o3a,  jp_fabm_o3ph, &
-                      jp_fabm_o3pc
+   INTEGER, PUBLIC :: jp_fabm_chl1,  jp_fabm_chl2, &
+                      jp_fabm_chl3,  jp_fabm_chl4, &
+                      jp_fabm_p1c,   jp_fabm_p1n,  &
+                      jp_fabm_p1p,   jp_fabm_p1s,  &
+                      jp_fabm_p2c,   jp_fabm_p2n,  &
+                      jp_fabm_p2p,   jp_fabm_p3c,  &
+                      jp_fabm_p3n,   jp_fabm_p3p,  &
+                      jp_fabm_p4c,   jp_fabm_p4n,  &
+                      jp_fabm_p4p,   jp_fabm_z4c,  &
+                      jp_fabm_z5c,   jp_fabm_z5n,  &
+                      jp_fabm_z5p,   jp_fabm_z6c,  &
+                      jp_fabm_z6n,   jp_fabm_z6p,  &
+                      jp_fabm_n1p,   jp_fabm_n3n,  &
+                      jp_fabm_n4n,   jp_fabm_n5s,  &
+                      jp_fabm_o2o,   jp_fabm_o3c,  &
+                      jp_fabm_o3ta,  jp_fabm_o3ba, &
+                      jp_fabm_o3pc,  jp_fabm_o3ph, &
+                      jp_fabm_r4n,   jp_fabm_r4c,  &
+                      jp_fabm_r4p,   jp_fabm_r6n,  &
+                      jp_fabm_r6c,   jp_fabm_r6p,  &
+                      jp_fabm_r6s,   jp_fabm_r8n,  &
+                      jp_fabm_r8c,   jp_fabm_r8p,  &
+                      jp_fabm_r8s,                 &
+                      jp_fabm_pgrow, jp_fabm_ploss
 
 #if defined key_fabm
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 10574)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90	(revision 10622)
@@ -105,9 +105,28 @@
       jp_fabm_o2o  = fabm_state_index( 'O2_o' )
       jp_fabm_o3c  = fabm_state_index( 'O3_c' )
-      jp_fabm_o3a  = fabm_state_index( 'O3_bioalk' )
+      jp_fabm_o3ba = fabm_state_index( 'O3_bioalk' )
+      jp_fabm_r4n  = fabm_state_index( 'R4_n' )
+      jp_fabm_r4c  = fabm_state_index( 'R4_c' )
+      jp_fabm_r4p  = fabm_state_index( 'R4_p' )
+      jp_fabm_r6n  = fabm_state_index( 'R6_n' )
+      jp_fabm_r6c  = fabm_state_index( 'R6_c' )
+      jp_fabm_r6p  = fabm_state_index( 'R6_p' )
+      jp_fabm_r6s  = fabm_state_index( 'R6_s' )
+      jp_fabm_r8n  = fabm_state_index( 'R8_n' )
+      jp_fabm_r8c  = fabm_state_index( 'R8_c' )
+      jp_fabm_r8p  = fabm_state_index( 'R8_p' )
+      jp_fabm_r8s  = fabm_state_index( 'R8_s' )
 
       ! Get indexes for select diagnostic variables
-      jp_fabm_o3ph = fabm_diag_index( 'O3_pH' )
-      jp_fabm_o3pc = fabm_diag_index( 'O3_pCO2' )
+      jp_fabm_o3ta  = fabm_diag_index( 'O3_TA' )
+      jp_fabm_o3ph  = fabm_diag_index( 'O3_pH' )
+      jp_fabm_o3pc  = fabm_diag_index( 'O3_pCO2' )
+      jp_fabm_pgrow = fabm_diag_index( 'p_grow_sum_result' )
+      jp_fabm_ploss = fabm_diag_index( 'p_loss_sum_result' )
+      
+      MLD_MAX(:,:)   = 0.0
+      PGROW_AVG(:,:) = 0.0
+      PLOSS_AVG(:,:) = 0.0
+      PHYT_AVG(:,:)  = 0.0
 
       IF (lwp) THEN
@@ -445,5 +464,5 @@
       END DO
       IF (fabm_state_index == -1) THEN
-         CALL ctl_stop( 'Could not find '//TRIM(state_name)//' state variable' )
+         CALL ctl_warn( 'Could not find '//TRIM(state_name)//' state variable' )
       ELSE
          IF (lwp) WRITE(numout,*) 'Index for '//TRIM(state_name)//' is: ', fabm_state_index
@@ -477,5 +496,5 @@
       END DO
       IF (fabm_diag_index == -1) THEN
-         CALL ctl_stop( 'Could not find '//TRIM(diag_name)//' diagnostic' )
+         CALL ctl_warn( 'Could not find '//TRIM(diag_name)//' diagnostic' )
       ELSE
          IF (lwp) WRITE(numout,*) 'Index for '//TRIM(diag_name)//' is: ', fabm_diag_index
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 10574)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90	(revision 10622)
@@ -33,4 +33,7 @@
    USE inputs_fabm
    USE vertical_movement_fabm
+   USE zdfmxl
+   USE asmbgc, ONLY: mld_choice_bgc
+   USE lbclnk
 
    !USE fldread         !  time interpolation
@@ -113,4 +116,6 @@
 
       CALL st2d_fabm_nxt( kt )
+      
+      CALL asmdiags_fabm( kt )
 
       IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrfabm )
@@ -130,4 +135,88 @@
 
    END SUBROUTINE trc_sms_fabm
+   
+   SUBROUTINE asmdiags_fabm( kt )
+      INTEGER, INTENT(IN) :: kt
+      INTEGER :: ji,jj,jk,jkmax
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d, zmld
+	   
+      IF (kt == nittrc000) THEN
+         MLD_MAX(:,:) = 0.0
+      ENDIF
+      PGROW_AVG(:,:) = 0.0
+      PLOSS_AVG(:,:) = 0.0
+      PHYT_AVG(:,:)  = 0.0
+	     
+      pgrow_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_pgrow)
+      ploss_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_ploss)
+      
+      SELECT CASE( mld_choice_bgc )
+      CASE ( 1 )                   ! Turbocline/mixing depth [W points]
+         zmld(:,:) = hmld(:,:)
+      CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
+         zmld(:,:) = hmlp(:,:)
+      CASE ( 3 )                   ! Kara MLD [Interpolated]
+#if defined key_karaml
+         IF ( ln_kara ) THEN
+            zmld(:,:) = hmld_kara(:,:)
+         ELSE
+            CALL ctl_stop( ' Kara mixed layer requested for BGC assimilation,', &
+               &           ' but ln_kara=.false.' )
+         ENDIF
+#else
+         CALL ctl_stop( ' Kara mixed layer requested for BGC assimilation,', &
+            &           ' but is not defined' )
+#endif
+      CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
+         zmld(:,:) = hmld_tref(:,:)
+      CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
+         zmld(:,:) = hmlpt(:,:)
+      END SELECT
+	
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            !
+            jkmax = jpk-1
+            DO jk = jpk-1, 1, -1
+               IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                  & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                  zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                  jkmax = jk
+               ENDIF
+            END DO
+            !
+            DO jk = 1, jkmax
+               PHYT_AVG(ji,jj) = PHYT_AVG(ji,jj) + &
+                  &              trn(ji,jj,jk,jp_fabm_p1n) + &
+                  &              trn(ji,jj,jk,jp_fabm_p2n) + &
+                  &              trn(ji,jj,jk,jp_fabm_p3n) + &
+                  &              trn(ji,jj,jk,jp_fabm_p4n)
+               IF ( pgrow_3d(ji,jj,jk) .GT. 0.0 ) THEN
+                  PGROW_AVG(ji,jj) = PGROW_AVG(ji,jj) + &
+                     &               pgrow_3d(ji,jj,jk)
+               ENDIF
+               IF ( ploss_3d(ji,jj,jk) .GT. 0.0 ) THEN
+                  PLOSS_AVG(ji,jj) = PLOSS_AVG(ji,jj) + &
+                     &               ploss_3d(ji,jj,jk)
+               ENDIF
+            END DO
+           
+            PHYT_AVG(ji,jj)  = PHYT_AVG(ji,jj)  / REAL(jkmax)
+            PGROW_AVG(ji,jj) = PGROW_AVG(ji,jj) / REAL(jkmax)
+            PLOSS_AVG(ji,jj) = PLOSS_AVG(ji,jj) / REAL(jkmax)
+	
+            IF ( zmld(ji,jj) .GT. MLD_MAX(ji,jj) ) THEN
+               MLD_MAX(ji,jj) = zmld(ji,jj)
+            ENDIF
+            !
+         END DO
+      END DO
+      
+      PHYT_AVG(:,:)  = PHYT_AVG(:,:)  * tmask(:,:,1)
+      PGROW_AVG(:,:) = PGROW_AVG(:,:) * tmask(:,:,1)
+      PLOSS_AVG(:,:) = PLOSS_AVG(:,:) * tmask(:,:,1)
+      MLD_MAX(:,:)   = MLD_MAX(:,:)   * tmask(:,:,1)
+	   
+   END SUBROUTINE asmdiags_fabm
 
    SUBROUTINE compute_fabm()
Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 10574)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 10622)
@@ -225,4 +225,11 @@
 #endif
 
+#if defined key_fabm
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  PGROW_AVG  !: Phytoplankton growth for use in ASM code
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  PLOSS_AVG  !: Phytoplankton loss   for use in ASM code
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  PHYT_AVG   !: Phytoplankton        for use in ASM code
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  MLD_MAX    !: Maximum MLD          for use in ASM code
+#endif
+
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3.1 , NEMO Consortium (2010)
@@ -253,4 +260,6 @@
 ! FABM <<<+++
          &      ln_trc_sbc(jptra)     , ln_trc_cbc(jptra)     , ln_trc_obc(jptra)     ,       &
+         &      PGROW_AVG(jpi,jpj)    , PLOSS_AVG(jpi,jpj)    , PHYT_AVG(jpi,jpj)     ,       &
+         &      MLD_MAX(jpi,jpj)      ,                                                       &
 #endif
 #if defined key_bdy
