Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9295)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9296)
@@ -1294,8 +1294,8 @@
     ln_asmiau = .false.    !  Logical switch for Incremental Analysis Updating (IAU)
     ln_seaiceinc = .false. !  Logical switch for applying sea ice increments
-    ln_logchlinc = .false. !  Logical switch for applying logchl increments
-    ln_logchlbal = .false. !  Logical switch for applying logchl multivariate balancing
-    ln_fco2inc = .false.   !  Logical switch for applying fCO2 increments
-    ln_pco2inc = .false.   !  Logical switch for applying pCO2 increments
+    ln_slchltotinc = .false. !  Logical switch for applying slchltot increments
+    ln_slchltotbal = .false. !  Logical switch for applying slchltot multivariate balancing
+    ln_sfco2inc = .false.   !  Logical switch for applying sfCO2 increments
+    ln_spco2inc = .false.   !  Logical switch for applying spCO2 increments
     ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing
     nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1]
Index: anches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 9295)
+++ 	(revision )
@@ -1,134 +1,0 @@
-MODULE asmbal
-   !!======================================================================
-   !!                       ***  MODULE asmbal  ***
-   !! Assimilation balancing interface: Write to file the balancing increments
-   !!                                   calculated for biogeochemistry
-   !!======================================================================
-   !!----------------------------------------------------------------------
-   !!   'key_asminc' : Switch on the assimilation increment interface
-   !!----------------------------------------------------------------------
-   !!   asm_bal_wri  : Write out the background state
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE par_kind, ONLY: &   ! Precision variables
-      & wp
-   USE iom                 ! I/O module
-   USE asminc              ! Main assimilation increments module
-   USE asmpar              ! Parameters for the assimilation interface
-#if defined key_medusa && defined key_foam_medusa
-   USE par_medusa
-#elif defined key_hadocc
-   USE par_hadocc
-#endif
-
-   IMPLICIT NONE
-
-   !! * Routine accessibility
-   PRIVATE
-   PUBLIC asm_bal_wri   !: Write out the background state
-
-CONTAINS
-
-   SUBROUTINE asm_bal_wri( kt )
-      !!-----------------------------------------------------------------------
-      !!
-      !!                  ***  ROUTINE asm_bal_wri ***
-      !!
-      !! ** Purpose : Write to file the balancing increments
-      !!              calculated for biogeochemistry
-      !!
-      !! ** Method  : Write to file the balancing increments
-      !!              calculated for biogeochemistry
-      !!
-      !! ** Action  :
-      !!                   
-      !! References :
-      !!
-      !! History    :
-      !!        ! 2014-08 (D. Ford)  Initial version, based on asm_bkg_wri
-      !!-----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( IN ) :: kt        ! Current time-step
-      !! * Local declarations
-      CHARACTER(LEN=50) :: cl_asmbal     ! Filename (with extension)
-      LOGICAL           :: llok          ! Check if file exists
-      INTEGER           :: inum          ! File unit number
-      REAL(wp)          :: zdate         ! Date
-      !!-----------------------------------------------------------------------
-     
-      ! Set things up
-      zdate = REAL( ndastp )
-      WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
-
-      ! Check if file exists
-      INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
-      IF( .NOT. llok ) THEN
-         IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
-            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
-
-         ! Define the output file       
-         CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
-
-         ! Write the information
-         CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate   )
-
-         IF ( ln_logchlinc ) THEN
-#if defined key_medusa
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', logchl_balinc(:,:,:,jpchn) )
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', logchl_balinc(:,:,:,jpchd) )
-            IF ( ln_logchlbal ) THEN
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', logchl_balinc(:,:,:,jpphn) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', logchl_balinc(:,:,:,jpphd) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', logchl_balinc(:,:,:,jppds) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', logchl_balinc(:,:,:,jpzmi) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', logchl_balinc(:,:,:,jpzme) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', logchl_balinc(:,:,:,jpdin) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', logchl_balinc(:,:,:,jpsil) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', logchl_balinc(:,:,:,jpfer) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', logchl_balinc(:,:,:,jpdet) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', logchl_balinc(:,:,:,jpdtc) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', logchl_balinc(:,:,:,jpdic) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', logchl_balinc(:,:,:,jpalk) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', logchl_balinc(:,:,:,jpoxy) )
-            ENDIF
-#elif defined key_hadocc
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', logchl_balinc(:,:,:,jp_had_phy) )
-            IF ( ln_logchlbal ) THEN
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', logchl_balinc(:,:,:,jp_had_nut) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', logchl_balinc(:,:,:,jp_had_zoo) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', logchl_balinc(:,:,:,jp_had_pdn) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', logchl_balinc(:,:,:,jp_had_dic) )
-               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', logchl_balinc(:,:,:,jp_had_alk) )
-            ENDIF
-#endif
-         ENDIF
-
-         IF ( ln_pco2inc ) THEN
-#if defined key_medusa
-            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jpdic) )
-            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jpalk) )
-#elif defined key_hadocc
-            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )
-            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )
-#endif
-         ELSE IF ( ln_fco2inc ) THEN
-#if defined key_medusa
-            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jpdic) )
-            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jpalk) )
-#elif defined key_hadocc
-            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )
-            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )
-#endif
-         ENDIF
-
-         CALL iom_close( inum )
-      ELSE
-         CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
-            &           ' Therefore not writing out balancing increments at this timestep', &
-            &           ' Something has probably gone wrong somewhere' )
-         IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
-            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
-      ENDIF
-                                 
-   END SUBROUTINE asm_bal_wri
-END MODULE asmbal
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 9296)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 9296)
@@ -0,0 +1,916 @@
+MODULE asmbgc
+   !!======================================================================
+   !!                       ***  MODULE asmbgc  ***
+   !! Routines and declarations for biogeochemical assimilation
+   !!======================================================================
+   !! History :  3.6  ! 2018-02 (D. Ford)  Adapted from existing modules
+   !!----------------------------------------------------------------------
+   !! 'key_asminc'          : assimilation increment interface
+   !! 'key_top'             : passive tracers
+   !! 'key_hadocc'          : HadOCC model
+   !! 'key_medusa'          : MEDUSA model
+   !! 'key_foam_medusa'     : MEDUSA extras for FOAM OBS and ASM
+   !! 'key_roam'            : MEDUSA extras for carbonate cycle
+   !! 'key_karaml'          : Kara mixed layer depth
+   !!----------------------------------------------------------------------
+   !! asm_bgc_check_options : check biogeochemical assimilation options
+   !! asm_bgc_init_incs     : initialise bgc increments
+   !! asm_bgc_init_bkg      : initialise bgc background
+   !! asm_bgc_bal_wri       : write out balancing increments
+   !! asm_bgc_bkg_wri       : write out bgc background
+   !! slchltot_asm_inc      : apply the slchltot increment
+   !! spco2_asm_inc         : apply the pCO2/fCO2 increment
+   !!----------------------------------------------------------------------
+   USE par_kind,      ONLY: wp             ! kind parameters
+   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
+   USE iom                                 ! i/o
+   USE zdfmxl, ONLY :  &  
+!   &  hmld_tref,       &   
+#if defined key_karaml
+   &  hmld_kara,       &
+   &  ln_kara,         &
+#endif   
+   &  hmld,            & 
+   &  hmlp,            &
+   &  hmlpt 
+#if defined key_top
+   USE trc, ONLY: &
+      & trn,      &
+      & trb
+   USE par_trc, ONLY: &
+      & jptra
+#endif
+#if defined key_medusa && defined key_foam_medusa
+   USE asmlogchlbal_medusa, ONLY: &
+      & asm_logchl_bal_medusa
+   USE asmpco2bal, ONLY: &
+      & asm_pco2_bal
+   USE par_medusa
+   USE mocsy_mainmod     ! f2pCO2
+   USE sms_medusa, ONLY: pgrow_avg, &
+      &                  ploss_avg, &
+      &                  phyt_avg,  &
+      &                  mld_max
+#elif defined key_hadocc
+   USE asmlogchlbal_hadocc, ONLY: &
+      & asm_logchl_bal_hadocc
+   USE asmpco2bal, ONLY: &
+      & asm_pco2_bal
+   USE par_hadocc
+   USE trc, ONLY: pgrow_avg, &
+       &     ploss_avg, &
+       &     phyt_avg,  &
+       &     mld_max,   &
+       &     HADOCC_CHL
+   USE had_bgc_const, ONLY: cchl_p
+#endif
+   USE asmpar    ! c_asmbkg, c_asmbal, nitbkg_r, nitdin_r, nitiaustr_r, nitiaufin_r
+   USE oce              ! tsn
+
+   IMPLICIT NONE
+   PRIVATE                   
+
+   PUBLIC asm_bgc_check_options
+   PUBLIC asm_bgc_init_incs
+   PUBLIC asm_bgc_init_bkg
+   PUBLIC asm_bgc_bal_wri
+   PUBLIC asm_bgc_bkg_wri
+   PUBLIC slchltot_asm_inc
+   PUBLIC spco2_asm_inc
+
+   LOGICAL, PUBLIC :: ln_balwri = .FALSE.        !: No output of balancing incs
+   LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE.   !: No slchltot increment
+   LOGICAL, PUBLIC :: ln_slchltotbal = .FALSE.   !: No slchltot balancing
+   LOGICAL, PUBLIC :: ln_spco2inc = .FALSE.      !: No pCO2 increment
+   LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE.      !: No fCO2 increment
+
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc !: Increment to background slchltot
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc   !: Increment to background pCO2
+#if defined key_top
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: slchltot_balinc  !: Increment to BGC variables from logchl assim
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: spco2_balinc    !: Increment to BGC variables from pCO2 assim
+#endif
+#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg  !: Background phyto growth for logchl balancing
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg  !: Background phyto loss for logchl balancing
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg   !: Background phyto for logchl balancing
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg    !: Background max MLD for logchl balancing
+   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg     !: Background tracer state variables
+#endif
+#if defined key_hadocc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: chl_bkg        !: Background surface chlorophyll
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: cchl_p_bkg     !: Background surface carbon:chlorophyll
+#endif
+   REAL(wp), PUBLIC :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation
+                                      !: <= 0 implies no maximum applied (switch turned off)
+                                      !:  > 0 implies maximum absolute chl increment capped at this value
+
+   INTEGER, PUBLIC :: mld_choice_bgc    = 1   !: choice of mld criteria to use for biogeochemistry assimilation
+                                      !: 1) hmld      - Turbocline/mixing depth                           [W points]
+                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points]
+                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated]
+                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points]
+                                      !: 5) hmlpt     - Density criterion (0.01 kg/m^3 change from 10m)   [T points]
+
+CONTAINS
+
+   SUBROUTINE asm_bgc_check_options
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_check_options  ***
+      !!
+      !! ** Purpose :   check validity of biogeochemical assimilation options
+      !!
+      !! ** Method  :   check settings of logicals
+      !!
+      !! ** Action  :   call ctl_stop/ctl_warn if required
+      !!
+      !! References :   asm_inc_init
+      !!---------------------------------------------------------------------------
+
+      IF ( ( ln_balwri ).AND.( .NOT. ln_slchltotinc ).AND.( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN
+         CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', &
+            &           ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.')
+         ln_balwri = .FALSE.
+      ENDIF
+
+      IF ( ( ln_slchltotbal ).AND.( .NOT. ln_slchltotinc ) ) THEN
+         CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', &
+            &           ' Not assimilating logchl, so ln_slchltotbal will be set to .false.')
+         ln_slchltotbal = .FALSE.
+      ENDIF
+
+      IF ( ( ln_spco2inc ).AND.( ln_sfco2inc ) ) THEN
+         CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
+      ENDIF
+
+   END SUBROUTINE asm_bgc_check_options
+
+   !!======================================================================
+   !!======================================================================
+   !!======================================================================
+
+   SUBROUTINE asm_bgc_init_incs( knum )
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_init_incs  ***
+      !!
+      !! ** Purpose :   initialise bgc increments
+      !!
+      !! ** Method  :   allocate arrays and read increments from file
+      !!
+      !! ** Action  :   allocate arrays and read increments from file
+      !!
+      !! References :   asm_inc_init
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file
+      !!
+      !!---------------------------------------------------------------------------
+
+      IF ( ln_slchltotinc ) THEN
+         ALLOCATE( slchltot_bkginc(jpi,jpj) )
+         slchltot_bkginc(:,:) = 0.0
+#if defined key_top
+         ALLOCATE( slchltot_balinc(jpi,jpj,jpk,jptra) )
+         slchltot_balinc(:,:,:,:) = 0.0
+#endif
+      ENDIF
+      IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
+         ALLOCATE( spco2_bkginc(jpi,jpj) )
+         spco2_bkginc(:,:) = 0.0
+#if defined key_top
+         ALLOCATE( spco2_balinc(jpi,jpj,jpk,jptra) )
+         spco2_balinc(:,:,:,:) = 0.0
+#endif
+      ENDIF
+
+      IF ( ln_slchltotinc ) THEN
+         CALL iom_get( knum, jpdom_autoglo, 'bckinslchltot', slchltot_bkginc(:,:), 1 )
+         ! Apply the masks
+         slchltot_bkginc(:,:) = slchltot_bkginc(:,:) * tmask(:,:,1)
+         ! Set missing increments to 0.0 rather than 1e+20
+         ! to allow for differences in masks
+         WHERE( ABS( slchltot_bkginc(:,:) ) > 1.0e+10 ) slchltot_bkginc(:,:) = 0.0
+      ENDIF
+
+      IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
+         IF ( ln_spco2inc ) THEN
+            CALL iom_get( knum, jpdom_autoglo, 'bckinspco2', spco2_bkginc(:,:), 1 )
+         ELSE IF ( ln_sfco2inc ) THEN
+            CALL iom_get( knum, jpdom_autoglo, 'bckinsfco2', spco2_bkginc(:,:), 1 )
+         ENDIF
+         ! Apply the masks
+         spco2_bkginc(:,:) = spco2_bkginc(:,:) * tmask(:,:,1)
+         ! Set missing increments to 0.0 rather than 1e+20
+         ! to allow for differences in masks
+         WHERE( ABS( spco2_bkginc(:,:) ) > 1.0e+10 ) spco2_bkginc(:,:) = 0.0
+      ENDIF
+
+   END SUBROUTINE asm_bgc_init_incs
+
+   !!======================================================================
+   !!======================================================================
+   !!======================================================================
+
+   SUBROUTINE asm_bgc_init_bkg
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_init_bkg  ***
+      !!
+      !! ** Purpose :   initialise bgc background
+      !!
+      !! ** Method  :   allocate arrays and read background from file
+      !!
+      !! ** Action  :   allocate arrays and read background from file
+      !!
+      !! References :   asm_inc_init
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER :: inum   ! i/o unit of background file
+      INTEGER :: jt     ! loop counter
+      !!
+      !!---------------------------------------------------------------------------
+
+#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
+      IF ( ln_slchltotinc ) THEN
+
+         ALLOCATE( pgrow_avg_bkg(jpi,jpj)        )
+         ALLOCATE( ploss_avg_bkg(jpi,jpj)        )
+         ALLOCATE( phyt_avg_bkg(jpi,jpj)         )
+         ALLOCATE( mld_max_bkg(jpi,jpj)          )
+         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
+         pgrow_avg_bkg(:,:)  = 0.0
+         ploss_avg_bkg(:,:)  = 0.0
+         phyt_avg_bkg(:,:)   = 0.0
+         mld_max_bkg(:,:)    = 0.0
+         tracer_bkg(:,:,:,:) = 0.0
+
+#if defined key_hadocc
+         ALLOCATE( chl_bkg(jpi,jpj)    )
+         ALLOCATE( cchl_p_bkg(jpi,jpj) )
+         chl_bkg(:,:)    = 0.0
+         cchl_p_bkg(:,:) = 0.0
+#endif
+         
+         !--------------------------------------------------------------------
+         ! Read background variables for slchltot assimilation
+         ! Some only required if performing balancing
+         !--------------------------------------------------------------------
+
+         CALL iom_open( c_asmbkg, inum )
+
+#if defined key_hadocc
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl',  chl_bkg    )
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
+         chl_bkg(:,:)    = chl_bkg(:,:)    * tmask(:,:,1)
+         cchl_p_bkg(:,:) = cchl_p_bkg(:,:) * tmask(:,:,1)
+#elif defined key_medusa
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
+#endif
+         
+         IF ( ln_slchltotbal ) THEN
+
+            CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
+            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
+            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  )
+            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
+            pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
+            ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
+            phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1)
+            mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1)
+
+#if defined key_hadocc
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
+#elif defined key_medusa
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
+#endif
+         ELSE IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
+#if defined key_hadocc
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
+#elif defined key_medusa
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+#endif
+            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
+            mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
+         ENDIF
+
+         CALL iom_close( inum )
+         
+         DO jt = 1, jptra
+            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
+         END DO
+      
+      ELSE IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
+
+         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
+         ALLOCATE( mld_max_bkg(jpi,jpj)          )
+         tracer_bkg(:,:,:,:) = 0.0
+         mld_max_bkg(:,:)    = 0.0
+
+         CALL iom_open( c_asmbkg, inum )
+         
+#if defined key_hadocc
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
+#elif defined key_medusa
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+#endif
+         CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
+
+         CALL iom_close( inum )
+         
+         DO jt = 1, jptra
+            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
+         END DO
+         mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
+      
+      ENDIF
+#endif
+
+   END SUBROUTINE asm_bgc_init_bkg
+
+   !!======================================================================
+   !!======================================================================
+   !!======================================================================
+
+   SUBROUTINE asm_bgc_bal_wri( kt )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                  ***  ROUTINE asm_bgc_bal_wri ***
+      !!
+      !! ** Purpose : Write to file the balancing increments
+      !!              calculated for biogeochemistry
+      !!
+      !! ** Method  : Write to file the balancing increments
+      !!              calculated for biogeochemistry
+      !!
+      !! ** Action  :
+      !!                   
+      !! References :
+      !!
+      !! History    :
+      !!        ! 2014-08 (D. Ford)  Initial version, based on asm_bkg_wri
+      !!-----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( IN ) :: kt        ! Current time-step
+      !! * Local declarations
+      CHARACTER(LEN=50) :: cl_asmbal     ! Filename (with extension)
+      LOGICAL           :: llok          ! Check if file exists
+      INTEGER           :: inum          ! File unit number
+      REAL(wp)          :: zdate         ! Date
+      !!-----------------------------------------------------------------------
+     
+      ! Set things up
+      zdate = REAL( ndastp )
+      WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
+
+      ! Check if file exists
+      INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
+      IF( .NOT. llok ) THEN
+         IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
+            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
+
+         ! Define the output file       
+         CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
+
+         ! Write the information
+         CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate   )
+
+         IF ( ln_slchltotinc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', slchltot_balinc(:,:,:,jpchn) )
+            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', slchltot_balinc(:,:,:,jpchd) )
+            IF ( ln_slchltotbal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', slchltot_balinc(:,:,:,jpphn) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', slchltot_balinc(:,:,:,jpphd) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', slchltot_balinc(:,:,:,jppds) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', slchltot_balinc(:,:,:,jpzmi) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', slchltot_balinc(:,:,:,jpzme) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', slchltot_balinc(:,:,:,jpdin) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', slchltot_balinc(:,:,:,jpsil) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', slchltot_balinc(:,:,:,jpfer) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', slchltot_balinc(:,:,:,jpdet) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', slchltot_balinc(:,:,:,jpdtc) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', slchltot_balinc(:,:,:,jpdic) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', slchltot_balinc(:,:,:,jpalk) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', slchltot_balinc(:,:,:,jpoxy) )
+            ENDIF
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', slchltot_balinc(:,:,:,jp_had_phy) )
+            IF ( ln_slchltotbal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', slchltot_balinc(:,:,:,jp_had_nut) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', slchltot_balinc(:,:,:,jp_had_zoo) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', slchltot_balinc(:,:,:,jp_had_pdn) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', slchltot_balinc(:,:,:,jp_had_dic) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', slchltot_balinc(:,:,:,jp_had_alk) )
+            ENDIF
+#endif
+         ENDIF
+
+         IF ( ln_spco2inc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', spco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', spco2_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', spco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', spco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ELSE IF ( ln_sfco2inc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', spco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', spco2_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', spco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', spco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ENDIF
+
+         CALL iom_close( inum )
+      ELSE
+         CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
+            &           ' Therefore not writing out balancing increments at this timestep', &
+            &           ' Something has probably gone wrong somewhere' )
+         IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
+            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
+      ENDIF
+                                 
+   END SUBROUTINE asm_bgc_bal_wri
+
+   !!======================================================================
+   !!======================================================================
+   !!======================================================================
+
+   SUBROUTINE asm_bgc_bkg_wri( kt, knum )
+      !!---------------------------------------------------------------------------
+      !!                    ***  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
+      !!
+      !!---------------------------------------------------------------------------
+
+#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(:,:,1)     )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,1)         )
+#elif defined key_medusa && defined key_foam_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) )
+#endif
+
+   END SUBROUTINE asm_bgc_bkg_wri
+
+   !!======================================================================
+   !!======================================================================
+   !!======================================================================
+
+   SUBROUTINE slchltot_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE logchl_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the chlorophyll assimilation increments.
+      !!
+      !! ** Method  : Calculate increments to state variables using nitrogen
+      !!              balancing.
+      !!              Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!----------------------------------------------------------------------
+      INTEGER,  INTENT(IN) :: kt        ! Current time step
+      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
+      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
+      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
+      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
+      !
+      INTEGER  :: jk              ! Loop counter
+      INTEGER  :: it              ! Index
+      REAL(wp) :: zincwgt         ! IAU weight for current time step
+      REAL(wp) :: zincper         ! IAU interval in seconds
+      !!----------------------------------------------------------------------
+
+      IF ( kt <= nit000 ) THEN
+
+         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
+
+#if defined key_medusa && defined key_foam_medusa
+         CALL asm_logchl_bal_medusa( slchltot_bkginc, zincper, mld_choice_bgc, &
+            &                        rn_maxchlinc, ln_slchltotbal, ll_asmdin,  &
+            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
+            &                        phyt_avg_bkg, mld_max_bkg,              &
+            &                        tracer_bkg, slchltot_balinc )
+#elif defined key_hadocc
+         CALL asm_logchl_bal_hadocc( slchltot_bkginc, zincper, mld_choice_bgc, &
+            &                        rn_maxchlinc, ln_slchltotbal, ll_asmdin,  &
+            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
+            &                        phyt_avg_bkg, mld_max_bkg,              &
+            &                        chl_bkg, cchl_p_bkg,                    &
+            &                        tracer_bkg, slchltot_balinc )
+#else
+         CALL ctl_stop( 'Attempting to assimilate slchltot, ', &
+            &           'but not defined a biogeochemical model' )
+#endif
+
+      ENDIF
+
+      IF ( ll_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
+            ! note this is not a tendency so should not be divided by rdt
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
+                  &  kt,' with IAU weight = ', pwgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the biogeochemical variables
+            ! Add directly to trn and trb, rather than to tra, because tra gets
+            ! reset to zero at the start of trc_stp, called after this routine
+#if defined key_medusa && defined key_foam_medusa
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          slchltot_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          slchltot_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+            END DO
+#elif defined key_hadocc
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
+                  &                          slchltot_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
+                  &                          slchltot_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+            END DO
+#endif
+
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+         ENDIF
+
+      ELSEIF ( ll_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+#if defined key_medusa && defined key_foam_medusa
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+               &                         slchltot_balinc(:,:,:,jp_msa0:jp_msa1)
+            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+#elif defined key_hadocc
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+               &                         slchltot_balinc(:,:,:,jp_had0:jp_had1)
+            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE slchltot_asm_inc
+
+   !!======================================================================
+   !!======================================================================
+   !!======================================================================
+
+   SUBROUTINE spco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
+      &                      ll_trainc, pt_bkginc, ps_bkginc )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE pco2_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
+      !!
+      !! ** Method  : Calculate increments to state variables using carbon
+      !!              balancing.
+      !!              Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(IN)                   :: kt           ! Current time step
+      LOGICAL, INTENT(IN)                   :: ll_asmdin    ! Flag for direct initialisation
+      LOGICAL, INTENT(IN)                   :: ll_asmiau    ! Flag for incremental analysis update
+      INTEGER, INTENT(IN)                   :: kcycper      ! Dimension of pwgtiau
+      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
+      LOGICAL, INTENT(IN)                   :: ll_trainc    ! Flag for T/S increments
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
+      !
+      INTEGER                               :: ji, jj, jk   ! Loop counters
+      INTEGER                               :: it           ! Index
+      INTEGER                               :: jkmax        ! Index of mixed layer depth
+      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
+      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
+      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
+      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
+      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
+      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
+      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
+      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
+      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
+
+      ! Coefficients for fCO2 to pCO2 conversion
+      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
+      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
+      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
+      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
+      REAL(wp)                              :: zcoef_fco2_5 = 2.0
+      REAL(wp)                              :: zcoef_fco2_6 = 57.7
+      REAL(wp)                              :: zcoef_fco2_7 = 0.118
+      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
+      !!----------------------------------------------------------------------
+
+      IF ( kt <= nit000 ) THEN
+
+         !--------------------------------------------------------------------
+         ! DIC and alkalinity balancing
+         !--------------------------------------------------------------------
+
+         IF ( ln_sfco2inc ) THEN
+#if defined key_medusa && defined key_roam
+            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
+            patm(1) = 1.0
+            phyd(1) = 0.0
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  CALL f2pCO2( spco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
+               END DO
+            END DO
+#else
+            ! If assimilating fCO2, then convert to pCO2 using temperature
+            ! See flux_gas.F90 within HadOCC for details of calculation
+            pco2_bkginc_temp(:,:) = spco2_bkginc(:,:) /                                                             &
+               &                    EXP((zcoef_fco2_1                                                            + &
+               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
+               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
+               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
+               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
+               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
+#endif
+         ELSE
+            pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
+         ENDIF
+
+         ! Account for physics assimilation if required
+         IF ( ll_trainc ) THEN
+            tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
+            sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
+         ELSE
+            tem_bkg_temp(:,:) = tsn(:,:,1,1)
+            sal_bkg_temp(:,:) = tsn(:,:,1,2)
+         ENDIF
+
+#if defined key_medusa
+         ! Account for logchl balancing if required
+         IF ( ln_slchltotinc .AND. ln_slchltotbal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + slchltot_balinc(:,:,1,jpdic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + slchltot_balinc(:,:,1,jpalk)
+         ELSE
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
+         ENDIF
+
+         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
+            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
+            &               spco2_balinc(:,:,1,jpdic), spco2_balinc(:,:,1,jpalk) )
+
+#elif defined key_hadocc
+         ! Account for slchltot balancing if required
+         IF ( ln_slchltotinc .AND. ln_slchltotbal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + slchltot_balinc(:,:,1,jp_had_dic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + slchltot_balinc(:,:,1,jp_had_alk)
+         ELSE
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
+         ENDIF
+
+         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
+            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
+            &               spco2_balinc(:,:,1,jp_had_dic), spco2_balinc(:,:,1,jp_had_alk) )
+
+#else
+         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
+            &           'but not defined a biogeochemical model' )
+#endif
+
+         ! Select mixed layer
+         IF ( ll_asmdin ) THEN
+#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
+            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
+               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
+            zmld(:,:) = mld_max_bkg(:,:)
+#else
+            CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
+#endif
+         ELSE
+            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 pCO2 assimilation,', &
+                     &           ' but ln_kara=.false.' )
+               ENDIF
+#else
+               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
+                  &           ' but is not defined' )
+#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' )
+            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
+               zmld(:,:) = hmlpt(:,:)
+            END SELECT
+         ENDIF
+
+         ! Propagate through mixed layer
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               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
+               !
+#if defined key_top
+               DO jk = 2, jkmax
+                  spco2_balinc(ji,jj,jk,:) = spco2_balinc(ji,jj,1,:)
+               END DO
+#endif
+               !
+            END DO
+         END DO
+
+      ENDIF
+
+      IF ( ll_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
+            ! note this is not a tendency so should not be divided by rdt
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               IF ( ln_spco2inc ) THEN
+                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
+                     &  kt,' with IAU weight = ', pwgtiau(it)
+               ELSE IF ( ln_sfco2inc ) THEN
+                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
+                     &  kt,' with IAU weight = ', pwgtiau(it)
+               ENDIF
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the biogeochemical variables
+            ! Add directly to trn and trb, rather than to tra, because tra gets
+            ! reset to zero at the start of trc_stp, called after this routine
+#if defined key_medusa && defined key_foam_medusa
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          spco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          spco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+            END DO
+#elif defined key_hadocc
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
+                  &                          spco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
+                  &                          spco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+            END DO
+#endif
+
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+
+         ENDIF
+
+      ELSEIF ( ll_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+#if defined key_medusa && defined key_foam_medusa
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+               &                         spco2_balinc(:,:,:,jp_msa0:jp_msa1)
+            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+#elif defined key_hadocc
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+               &                         spco2_balinc(:,:,:,jp_had0:jp_had1)
+            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE spco2_asm_inc
+
+   !!======================================================================
+
+END MODULE asmbgc
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 9295)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 9296)
@@ -50,20 +50,6 @@
    USE ice
 #endif
-#if defined key_hadocc
-   USE trc, ONLY: trn, &
-      &     pgrow_avg, &
-      &     ploss_avg, &
-      &     phyt_avg,  &
-      &     mld_max,   &
-      &     HADOCC_CHL
-   USE had_bgc_const, ONLY: cchl_p
-   USE par_hadocc
-#elif defined key_medusa && defined key_foam_medusa
-   USE trc, ONLY: trn
-   USE sms_medusa, ONLY: pgrow_avg, &
-      &                  ploss_avg, &
-      &                  phyt_avg,  &
-      &                  mld_max
-   USE par_medusa
+#if defined key_top
+   USE asmbgc, ONLY: asm_bgc_bkg_wri
 #endif
    IMPLICIT NONE
@@ -138,39 +124,9 @@
 !            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               )
             CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
-#if defined key_hadocc
-            CALL iom_rstput( kt, nitbkg_r, inum, 'pgrow_avg'   , pgrow_avg             )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'ploss_avg'   , ploss_avg             )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'phyt_avg'    , phyt_avg              )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'mld_max'     , mld_max               )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_zoo'  , trn(:,:,:,jp_had_zoo) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_pdn'  , trn(:,:,:,jp_had_pdn) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_dic'  , trn(:,:,:,jp_had_dic) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_alk'  , trn(:,:,:,jp_had_alk) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_chl'  , HADOCC_CHL(:,:,1)     )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_cchl' , cchl_p(:,:,1)         )
-#elif defined key_medusa && defined key_foam_medusa
-            CALL iom_rstput( kt, nitbkg_r, inum, 'pgrow_avg'   , pgrow_avg        )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'ploss_avg'   , ploss_avg        )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'phyt_avg'    , phyt_avg         )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'mld_max'     , mld_max          )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_chn'  , trn(:,:,:,jpchn) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_chd'  , trn(:,:,:,jpchd) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_phn'  , trn(:,:,:,jpphn) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_phd'  , trn(:,:,:,jpphd) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_pds'  , trn(:,:,:,jppds) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_zmi'  , trn(:,:,:,jpzmi) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_zme'  , trn(:,:,:,jpzme) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_din'  , trn(:,:,:,jpdin) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_sil'  , trn(:,:,:,jpsil) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_fer'  , trn(:,:,:,jpfer) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_det'  , trn(:,:,:,jpdet) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_dtc'  , trn(:,:,:,jpdtc) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_dic'  , trn(:,:,:,jpdic) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_alk'  , trn(:,:,:,jpalk) )
-            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_oxy'  , trn(:,:,:,jpoxy) )
+            !
+#if defined key_top
+            CALL asm_bgc_bkg_wri( kt, inum )
+            !
 #endif
-            !
             CALL iom_close( inum )
          ENDIF
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 9295)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 9296)
@@ -22,6 +22,4 @@
    !!   ssh_asm_inc    : Apply the SSH increment
    !!   seaice_asm_inc : Apply the seaice increment
-   !!   logchl_asm_inc : Apply the logchl increment
-   !!   pco2_asm_inc   : Apply the pCO2/fCO2 increment
    !!----------------------------------------------------------------------
    USE wrk_nemo         ! Memory Allocation
@@ -46,34 +44,5 @@
 #endif
    USE sbc_oce          ! Surface boundary condition variables.
-   USE zdfmxl, ONLY :  &  
-!   &  hmld_tref,       &   
-#if defined key_karaml
-   &  hmld_kara,       &
-   &  ln_kara,         &
-#endif   
-   &  hmld,            & 
-   &  hmlp,            &
-   &  hmlpt 
-#if defined key_top
-   USE trc, ONLY: &
-      & trn,      &
-      & trb
-   USE par_trc, ONLY: &
-      & jptra
-#endif
-#if defined key_medusa && defined key_foam_medusa
-   USE asmlogchlbal_medusa, ONLY: &
-      & asm_logchl_bal_medusa
-   USE asmpco2bal, ONLY: &
-      & asm_pco2_bal
-   USE par_medusa
-   USE mocsy_mainmod
-#elif defined key_hadocc
-   USE asmlogchlbal_hadocc, ONLY: &
-      & asm_logchl_bal_hadocc
-   USE asmpco2bal, ONLY: &
-      & asm_pco2_bal
-   USE par_hadocc
-#endif
+   USE asmbgc           ! Biogeochemistry assimilation
 
    IMPLICIT NONE
@@ -86,6 +55,5 @@
    PUBLIC   ssh_asm_inc    !: Apply the SSH increment
    PUBLIC   seaice_asm_inc !: Apply the seaice increment
-   PUBLIC   logchl_asm_inc !: Apply the logchl increment
-   PUBLIC   pco2_asm_inc   !: Apply the pCO2/fCO2 increment
+   PUBLIC   bgc_asm_inc    !: Apply the biogeochemistry increments
 
 #if defined key_asminc
@@ -95,5 +63,4 @@
 #endif
    LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields
-   LOGICAL, PUBLIC :: ln_balwri = .FALSE.      !: No output of the assimilation balancing increments
    LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment
    LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization
@@ -102,8 +69,5 @@
    LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
    LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment
-   LOGICAL, PUBLIC :: ln_logchlinc = .FALSE.   !: No log10(chlorophyll) increment
-   LOGICAL, PUBLIC :: ln_logchlbal = .FALSE.   !: Don't apply multivariate log10(chlorophyll) balancing
-   LOGICAL, PUBLIC :: ln_pco2inc = .FALSE.     !: No pCO2 increment
-   LOGICAL, PUBLIC :: ln_fco2inc = .FALSE.     !: No fCO2 increment
+   LOGICAL, PUBLIC :: lk_bgcinc = .FALSE.      !: No biogeochemistry increments
    LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
@@ -130,32 +94,4 @@
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
-
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pco2_bkginc   !: Increment to background pCO2
-#if defined key_top
-   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim
-   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc    !: Increment to BGC variables from pCO2 assim
-#endif
-#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
-   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg  !: Background phyto growth for logchl balancing
-   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg  !: Background phyto loss for logchl balancing
-   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg   !: Background phyto for logchl balancing
-   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg    !: Background max MLD for logchl balancing
-   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg     !: Background tracer state variables
-#endif
-#if defined key_hadocc
-   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: chl_bkg        !: Background surface chlorophyll
-   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: cchl_p_bkg     !: Background surface carbon:chlorophyll
-#endif
-   REAL(wp) :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation
-                                      !: <= 0 implies no maximum applied (switch turned off)
-                                      !:  > 0 implies maximum absolute chl increment capped at this value
-
-   INTEGER :: mld_choice_bgc    = 1   !: choice of mld criteria to use for biogeochemistry assimilation
-                                      !: 1) hmld      - Turbocline/mixing depth                           [W points]
-                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points]
-                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated]
-                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points]
-                                      !: 5) hmlpt     - Density criterion (0.01 kg/m^3 change from 10m)   [T points]
 
    !! * Substitutions
@@ -202,6 +138,6 @@
       NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri,                           &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
-         &                 ln_logchlinc, ln_logchlbal,                     &
-         &                 ln_pco2inc, ln_fco2inc,                         &
+         &                 ln_slchltotinc, ln_slchltotbal,                 &
+         &                 ln_spco2inc, ln_sfco2inc,                       &
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
@@ -239,8 +175,8 @@
          WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
          WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
-         WRITE(numout,*) '      Logical switch for applying logchl increments         ln_logchlinc = ', ln_logchlinc
-         WRITE(numout,*) '      Logical switch for logchl multivariate balancing      ln_logchlbal = ', ln_logchlbal
-         WRITE(numout,*) '      Logical switch for applying pCO2 increments             ln_pco2inc = ', ln_pco2inc
-         WRITE(numout,*) '      Logical switch for applying fCO2 increments             ln_fco2inc = ', ln_fco2inc
+         WRITE(numout,*) '      Logical switch for applying slchltot increments     ln_slchltotinc = ', ln_slchltotinc
+         WRITE(numout,*) '      Logical switch for slchltot multivariate balancing  ln_slchltotbal = ', ln_slchltotbal
+         WRITE(numout,*) '      Logical switch for applying pCO2 increments            ln_spco2inc = ', ln_spco2inc
+         WRITE(numout,*) '      Logical switch for applying fCO2 increments            ln_sfco2inc = ', ln_sfco2inc
          WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
          WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
@@ -289,4 +225,8 @@
          WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date
       ENDIF
+      
+      IF ( ln_slchltotinc .OR. ln_spco2inc .OR. ln_sfco2inc ) THEN
+         lk_bgcinc = .TRUE.
+      ENDIF
 
       IF ( nacc /= 0 ) &
@@ -301,7 +241,7 @@
       IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
          & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. &
-         &        ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) )) &
+         &        ( lk_bgcinc ) )) &
          & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
-         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc is set to .true.', &
+         &                ' ln_(bgc-variable)inc is set to .true.', &
          &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
          &                ' Inconsistent options')
@@ -312,7 +252,7 @@
 
       IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
-         & .AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) )  &
+         & .AND.( .NOT. lk_bgcinc ) )  &
          & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
-         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc are set to .false. :', &
+         &                ' ln_(bgc-variable)inc are set to .false. :', &
          &                ' The assimilation increments are not applied')
 
@@ -334,20 +274,6 @@
          &                ' Background time step for Direct Initialization is outside', &
          &                ' the cycle interval')
-
-      IF ( ( ln_balwri ).AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) ) THEN
-         CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', &
-            &           ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.')
-         ln_balwri = .FALSE.
-      ENDIF
-
-      IF ( ( ln_logchlbal ).AND.( .NOT. ln_logchlinc ) ) THEN
-         CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', &
-            &           ' Not assimilating logchl, so ln_logchlbal will be set to .false.')
-         ln_logchlbal = .FALSE.
-      ENDIF
-
-      IF ( ( ln_pco2inc ).AND.( ln_fco2inc ) ) THEN
-         CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
-      ENDIF
+      
+      IF ( lk_bgcinc ) CALL asm_bgc_check_options
 
       IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
@@ -446,22 +372,6 @@
       ssh_iau(:,:)    = 0.0
 #endif
-      IF ( ln_logchlinc ) THEN
-         ALLOCATE( logchl_bkginc(jpi,jpj) )
-         logchl_bkginc(:,:) = 0.0
-#if defined key_top
-         ALLOCATE( logchl_balinc(jpi,jpj,jpk,jptra) )
-         logchl_balinc(:,:,:,:) = 0.0
-#endif
-      ENDIF
-      IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
-         ALLOCATE( pco2_bkginc(jpi,jpj) )
-         pco2_bkginc(:,:) = 0.0
-#if defined key_top
-         ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
-         pco2_balinc(:,:,:,:) = 0.0
-#endif
-      ENDIF
       IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
-         &  .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN
+         &  .OR.( lk_bgcinc ) ) THEN
 
          !--------------------------------------------------------------------
@@ -537,24 +447,6 @@
          ENDIF
 
-         IF ( ln_logchlinc ) THEN
-            CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl', logchl_bkginc(:,:), 1 )
-            ! Apply the masks
-            logchl_bkginc(:,:) = logchl_bkginc(:,:) * tmask(:,:,1)
-            ! Set missing increments to 0.0 rather than 1e+20
-            ! to allow for differences in masks
-            WHERE( ABS( logchl_bkginc(:,:) ) > 1.0e+10 ) logchl_bkginc(:,:) = 0.0
-         ENDIF
-
-         IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
-            IF ( ln_pco2inc ) THEN
-               CALL iom_get( inum, jpdom_autoglo, 'bckinpco2', pco2_bkginc(:,:), 1 )
-            ELSE IF ( ln_fco2inc ) THEN
-               CALL iom_get( inum, jpdom_autoglo, 'bckinfco2', pco2_bkginc(:,:), 1 )
-            ENDIF
-            ! Apply the masks
-            pco2_bkginc(:,:) = pco2_bkginc(:,:) * tmask(:,:,1)
-            ! Set missing increments to 0.0 rather than 1e+20
-            ! to allow for differences in masks
-            WHERE( ABS( pco2_bkginc(:,:) ) > 1.0e+10 ) pco2_bkginc(:,:) = 0.0
+         IF ( lk_bgcinc ) THEN
+            CALL asm_bgc_init_incs( inum )
          ENDIF
 
@@ -670,121 +562,8 @@
 
       ENDIF
-
-#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
-      IF ( ln_logchlinc ) THEN
-
-         ALLOCATE( pgrow_avg_bkg(jpi,jpj)        )
-         ALLOCATE( ploss_avg_bkg(jpi,jpj)        )
-         ALLOCATE( phyt_avg_bkg(jpi,jpj)         )
-         ALLOCATE( mld_max_bkg(jpi,jpj)          )
-         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
-         pgrow_avg_bkg(:,:)  = 0.0
-         ploss_avg_bkg(:,:)  = 0.0
-         phyt_avg_bkg(:,:)   = 0.0
-         mld_max_bkg(:,:)    = 0.0
-         tracer_bkg(:,:,:,:) = 0.0
-
-#if defined key_hadocc
-         ALLOCATE( chl_bkg(jpi,jpj)    )
-         ALLOCATE( cchl_p_bkg(jpi,jpj) )
-         chl_bkg(:,:)    = 0.0
-         cchl_p_bkg(:,:) = 0.0
-#endif
          
-         !--------------------------------------------------------------------
-         ! Read background variables for logchl assimilation
-         ! Some only required if performing balancing
-         !--------------------------------------------------------------------
-
-         CALL iom_open( c_asmbkg, inum )
-
-#if defined key_hadocc
-         CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl',  chl_bkg    )
-         CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
-         chl_bkg(:,:)    = chl_bkg(:,:)    * tmask(:,:,1)
-         cchl_p_bkg(:,:) = cchl_p_bkg(:,:) * tmask(:,:,1)
-#elif defined key_medusa
-         CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
-         CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
-#endif
-         
-         IF ( ln_logchlbal ) THEN
-
-            CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
-            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
-            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  )
-            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
-            pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
-            ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
-            phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1)
-            mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1)
-
-#if defined key_hadocc
-            CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
-            CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
-            CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
-            CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
-            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
-            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
-#elif defined key_medusa
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
-#endif
-         ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
-#if defined key_hadocc
-            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
-            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
-#elif defined key_medusa
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
-            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
-#endif
-            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
-            mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
-         ENDIF
-
-         CALL iom_close( inum )
-         
-         DO jt = 1, jptra
-            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
-         END DO
-      
-      ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
-
-         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
-         ALLOCATE( mld_max_bkg(jpi,jpj)          )
-         tracer_bkg(:,:,:,:) = 0.0
-         mld_max_bkg(:,:)    = 0.0
-
-         CALL iom_open( c_asmbkg, inum )
-         
-#if defined key_hadocc
-         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
-         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
-#elif defined key_medusa
-         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
-         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
-#endif
-         CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
-
-         CALL iom_close( inum )
-         
-         DO jt = 1, jptra
-            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
-         END DO
-         mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
-      
-      ENDIF
-#endif
+      IF ( lk_bgcinc ) THEN
+         CALL asm_bgc_init_bkg
+      ENDIF
       !
    END SUBROUTINE asm_inc_init
@@ -1408,381 +1187,35 @@
    END SUBROUTINE seaice_asm_inc
 
-   SUBROUTINE logchl_asm_inc( kt )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE logchl_asm_inc  ***
+
+   SUBROUTINE bgc_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE bgc_asm_inc  ***
       !!          
-      !! ** Purpose : Apply the chlorophyll assimilation increments.
-      !!
-      !! ** Method  : Calculate increments to state variables using nitrogen
-      !!              balancing.
-      !!              Direct initialization or Incremental Analysis Updating.
-      !!
-      !! ** Action  : 
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(IN) :: kt   ! Current time step
-      !
-      INTEGER  :: jk              ! Loop counter
-      INTEGER  :: it              ! Index
-      REAL(wp) :: zincwgt         ! IAU weight for current time step
-      REAL(wp) :: zincper         ! IAU interval in seconds
-      !!----------------------------------------------------------------------
-
-      IF ( kt <= nit000 ) THEN
-
-         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
-
-#if defined key_medusa && defined key_foam_medusa
-         CALL asm_logchl_bal_medusa( logchl_bkginc, zincper, mld_choice_bgc, &
-            &                        rn_maxchlinc, ln_logchlbal, ln_asmdin,  &
-            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
-            &                        phyt_avg_bkg, mld_max_bkg,              &
-            &                        tracer_bkg, logchl_balinc )
-#elif defined key_hadocc
-         CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, &
-            &                        rn_maxchlinc, ln_logchlbal, ln_asmdin,  &
-            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
-            &                        phyt_avg_bkg, mld_max_bkg,              &
-            &                        chl_bkg, cchl_p_bkg,                    &
-            &                        tracer_bkg, logchl_balinc )
-#else
-         CALL ctl_stop( 'Attempting to assimilate logchl, ', &
-            &           'but not defined a biogeochemical model' )
-#endif
-
-      ENDIF
-
-      IF ( ln_asmiau ) THEN
-
-         !--------------------------------------------------------------------
-         ! Incremental Analysis Updating
-         !--------------------------------------------------------------------
-
-         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
-
-            it = kt - nit000 + 1
-            zincwgt = wgtiau(it)      ! IAU weight for the current time step
-            ! note this is not a tendency so should not be divided by rdt
-
-            IF(lwp) THEN
-               WRITE(numout,*) 
-               WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
-                  &  kt,' with IAU weight = ', wgtiau(it)
-               WRITE(numout,*) '~~~~~~~~~~~~'
-            ENDIF
-
-            ! Update the biogeochemical variables
-            ! Add directly to trn and trb, rather than to tra, because tra gets
-            ! reset to zero at the start of trc_stp, called after this routine
-#if defined key_medusa && defined key_foam_medusa
-            DO jk = 1, jpkm1
-               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
-                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
-               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
-                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
-            END DO
-#elif defined key_hadocc
-            DO jk = 1, jpkm1
-               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
-                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
-               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
-                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
-            END DO
-#endif
-
-            ! Do not deallocate arrays - needed by asm_bal_wri
-            ! which is called at end of model run
-         ENDIF
-
-      ELSEIF ( ln_asmdin ) THEN 
-
-         !--------------------------------------------------------------------
-         ! Direct Initialization
-         !--------------------------------------------------------------------
-         
-         IF ( kt == nitdin_r ) THEN
-
-            neuler = 0                    ! Force Euler forward step
-
-#if defined key_medusa && defined key_foam_medusa
-            ! Initialize the now fields with the background + increment
-            ! Background currently is what the model is initialised with
-            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
-               &           ' Background state is taken from model rather than background file' )
-            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
-               &                         logchl_balinc(:,:,:,jp_msa0:jp_msa1)
-            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
-#elif defined key_hadocc
-            ! Initialize the now fields with the background + increment
-            ! Background currently is what the model is initialised with
-            CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
-               &           ' Background state is taken from model rather than background file' )
-            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
-               &                         logchl_balinc(:,:,:,jp_had0:jp_had1)
-            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
-#endif
- 
-            ! Do not deallocate arrays - needed by asm_bal_wri
-            ! which is called at end of model run
-         ENDIF
-         !
-      ENDIF
-      !
-   END SUBROUTINE logchl_asm_inc
-
-
-   SUBROUTINE pco2_asm_inc( kt )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE pco2_asm_inc  ***
-      !!          
-      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
-      !!
-      !! ** Method  : Calculate increments to state variables using carbon
-      !!              balancing.
-      !!              Direct initialization or Incremental Analysis Updating.
-      !!
-      !! ** Action  : 
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(IN)                   :: kt           ! Current time step
-      !
-      INTEGER                               :: ji, jj, jk   ! Loop counters
-      INTEGER                               :: it           ! Index
-      INTEGER                               :: jkmax        ! Index of mixed layer depth
-      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
-      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
-      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
-      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
-      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
-      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
-      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
-      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
-      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
-
-      ! Coefficients for fCO2 to pCO2 conversion
-      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
-      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
-      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
-      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
-      REAL(wp)                              :: zcoef_fco2_5 = 2.0
-      REAL(wp)                              :: zcoef_fco2_6 = 57.7
-      REAL(wp)                              :: zcoef_fco2_7 = 0.118
-      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
-      !!----------------------------------------------------------------------
-
-      IF ( kt <= nit000 ) THEN
-
-         !--------------------------------------------------------------------
-         ! DIC and alkalinity balancing
-         !--------------------------------------------------------------------
-
-         IF ( ln_fco2inc ) THEN
-#if defined key_medusa && defined key_roam
-            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
-            patm(1) = 1.0
-            phyd(1) = 0.0
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  CALL f2pCO2( pco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
-               END DO
-            END DO
-#else
-            ! If assimilating fCO2, then convert to pCO2 using temperature
-            ! See flux_gas.F90 within HadOCC for details of calculation
-            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) /                                                             &
-               &                    EXP((zcoef_fco2_1                                                            + &
-               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
-               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
-               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
-               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
-               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
-#endif
-         ELSE
-            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:)
-         ENDIF
-
-         ! Account for physics assimilation if required
-         IF ( ln_trainc ) THEN
-            tem_bkg_temp(:,:) = tsn(:,:,1,1) + t_bkginc(:,:,1)
-            sal_bkg_temp(:,:) = tsn(:,:,1,2) + s_bkginc(:,:,1)
-         ELSE
-            tem_bkg_temp(:,:) = tsn(:,:,1,1)
-            sal_bkg_temp(:,:) = tsn(:,:,1,2)
-         ENDIF
-
-#if defined key_medusa
-         ! Account for logchl balancing if required
-         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN
-            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + logchl_balinc(:,:,1,jpdic)
-            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + logchl_balinc(:,:,1,jpalk)
-         ELSE
-            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
-            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
-         ENDIF
-
-         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
-            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
-            &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
-
-#elif defined key_hadocc
-         ! Account for logchl balancing if required
-         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN
-            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + logchl_balinc(:,:,1,jp_had_dic)
-            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + logchl_balinc(:,:,1,jp_had_alk)
-         ELSE
-            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
-            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
-         ENDIF
-
-         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
-            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
-            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
-
-#else
-         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
-            &           'but not defined a biogeochemical model' )
-#endif
-
-         ! Select mixed layer
-         IF ( ln_asmdin ) THEN
-#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
-            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
-               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
-            zmld(:,:) = mld_max_bkg(:,:)
-#else
-            CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
-#endif
-         ELSE
-            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 pCO2 assimilation,', &
-                     &           ' but ln_kara=.false.' )
-               ENDIF
-#else
-               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
-                  &           ' but is not defined' )
-#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' )
-            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
-               zmld(:,:) = hmlpt(:,:)
-            END SELECT
-         ENDIF
-
-         ! Propagate through mixed layer
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               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
-               !
-#if defined key_top
-               DO jk = 2, jkmax
-                  pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
-               END DO
-#endif
-               !
-            END DO
-         END DO
-
-      ENDIF
-
-      IF ( ln_asmiau ) THEN
-
-         !--------------------------------------------------------------------
-         ! Incremental Analysis Updating
-         !--------------------------------------------------------------------
-
-         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
-
-            it = kt - nit000 + 1
-            zincwgt = wgtiau(it)      ! IAU weight for the current time step
-            ! note this is not a tendency so should not be divided by rdt
-
-            IF(lwp) THEN
-               WRITE(numout,*) 
-               IF ( ln_pco2inc ) THEN
-                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
-                     &  kt,' with IAU weight = ', wgtiau(it)
-               ELSE IF ( ln_fco2inc ) THEN
-                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
-                     &  kt,' with IAU weight = ', wgtiau(it)
-               ENDIF
-               WRITE(numout,*) '~~~~~~~~~~~~'
-            ENDIF
-
-            ! Update the biogeochemical variables
-            ! Add directly to trn and trb, rather than to tra, because tra gets
-            ! reset to zero at the start of trc_stp, called after this routine
-#if defined key_medusa && defined key_foam_medusa
-            DO jk = 1, jpkm1
-               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
-                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
-               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
-                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
-            END DO
-#elif defined key_hadocc
-            DO jk = 1, jpkm1
-               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
-                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
-               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
-                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
-            END DO
-#endif
-
-            ! Do not deallocate arrays - needed by asm_bal_wri
-            ! which is called at end of model run
-
-         ENDIF
-
-      ELSEIF ( ln_asmdin ) THEN 
-
-         !--------------------------------------------------------------------
-         ! Direct Initialization
-         !--------------------------------------------------------------------
-         
-         IF ( kt == nitdin_r ) THEN
-
-            neuler = 0                    ! Force Euler forward step
-
-#if defined key_medusa && defined key_foam_medusa
-            ! Initialize the now fields with the background + increment
-            ! Background currently is what the model is initialised with
-            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
-               &           ' Background state is taken from model rather than background file' )
-            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
-               &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1)
-            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
-#elif defined key_hadocc
-            ! Initialize the now fields with the background + increment
-            ! Background currently is what the model is initialised with
-            CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
-               &           ' Background state is taken from model rather than background file' )
-            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
-               &                         pco2_balinc(:,:,:,jp_had0:jp_had1)
-            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
-#endif
- 
-            ! Do not deallocate arrays - needed by asm_bal_wri
-            ! which is called at end of model run
-         ENDIF
-         !
-      ENDIF
-      !
-   END SUBROUTINE pco2_asm_inc
+      !! ** Purpose : Apply the biogeochemistry assimilation increments
+      !!
+      !! ** Method  : Call relevant routines in asmbgc
+      !!
+      !! ** Action  : Call relevant routines in asmbgc
+      !!
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER, INTENT(in   ) :: kt        ! Current time step
+      !
+      INTEGER                :: icycper   ! Dimension of wgtiau
+      !!
+      !!----------------------------------------------------------------------
+      
+      icycper = SIZE( wgtiau )
+      
+      IF( ln_slchltotinc ) THEN
+         CALL slchltot_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau )
+      ENDIF
+      
+      IF( ln_sfco2inc .OR. ln_spco2inc ) THEN
+         CALL spco2_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, &
+            &                ln_trainc, t_bkginc, s_bkginc )
+      ENDIF
+
+   END SUBROUTINE bgc_asm_inc
    
    !!======================================================================
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 9295)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 9296)
@@ -61,5 +61,5 @@
    USE asminc          ! assimilation increments     
    USE asmbkg          ! writing out state trajectory
-   USE asmbal          ! writing out biogeochemical assimilation balancing increments
+   USE asmbgc          ! biogeochemical assimilation increments
    USE diaptr          ! poleward transports           (dia_ptr_init routine)
    USE diadct          ! sections transports           (dia_dct_init routine)
@@ -161,6 +161,5 @@
                 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics
                 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
-                IF( ln_logchlinc ) CALL logchl_asm_inc( nit000 - 1 )
-                IF( ln_fco2inc .OR. ln_pco2inc ) CALL pco2_asm_inc( nit000 - 1 )
+                IF( lk_bgcinc ) CALL bgc_asm_inc( nit000 - 1 )    ! BGC
              ENDIF
           ENDIF
@@ -195,5 +194,5 @@
       IF( lk_diaobs   )   CALL dia_obs_wri
       !
-      IF( ( lk_asminc ).AND.( ln_balwri ) ) CALL asm_bal_wri( nitend )  ! Output balancing increments
+      IF( ( lk_asminc ).AND.( ln_balwri ) ) CALL asm_bgc_bal_wri( nitend )  ! Output balancing increments
       !
       IF( ln_icebergs )   CALL icb_end( nitend )
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 9295)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 9296)
@@ -253,8 +253,6 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_asminc .AND. ln_asmiau .AND. ln_logchlinc ) &
-         &               CALL logchl_asm_inc( kstp )  ! logchl assimilation
-      IF( lk_asminc .AND. ln_asmiau .AND. ( ln_fco2inc .OR. ln_pco2inc ) ) &
-         &               CALL pco2_asm_inc( kstp )    ! pCO2/fCO2 assimilation
+      IF( lk_asminc .AND. ln_asmiau .AND. lk_bgcinc ) &
+         &               CALL bgc_asm_inc( kstp )     ! biogeochemistry assimilation
                          CALL trc_stp( kstp )         ! time-stepping
 #endif
