Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_medusa_ref
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_medusa_ref	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_medusa_ref	(revision 10149)
@@ -305,4 +305,7 @@
       scl_chl    = 1.0      !! (IF key_oasis) scaling the exported to UM chlorophyll field
 !!
+!! FOAM                         !! DAF
+      ln_foam_medusa = .false.  !! Calculate and save diagnostics for data assimilation
+!!
 /
 !'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
@@ -319,4 +322,5 @@
 !!                    = (119 / 106) = 1.1226 mol O2 / mol C
 !!       xo2min    :  oxygen minimum concentration
+!!       xobs_xco2a:  observed atmospheric xCO2 (not used if <= 0)
 !
 &natroam
@@ -326,4 +330,5 @@
    xthetarem = 1.1226
    xo2min    = 4.0
+   xobs_xco2a= -999.0
 /
 !'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 10149)
@@ -1287,20 +1287,43 @@
 &nam_asminc   !   assimilation increments                               ('key_asminc')
 !-----------------------------------------------------------------------
-    ln_bkgwri = .false.    !  Logical switch for writing out background state
-    ln_trainc = .false.    !  Logical switch for applying tracer increments
-    ln_dyninc = .false.    !  Logical switch for applying velocity increments
-    ln_sshinc = .false.    !  Logical switch for applying SSH increments
-    ln_asmdin = .false.    !  Logical switch for Direct Initialization (DI)
-    ln_asmiau = .false.    !  Logical switch for Incremental Analysis Updating (IAU)
-    ln_seaiceinc = .false. !  Logical switch for applying sea ice increments
+    ln_bkgwri      = .false. !  Logical switch for writing out background state
+    ln_balwri      = .false. !  Logical switch for writing out balancing increments
+    ln_trainc      = .false. !  Logical switch for applying tracer increments
+    ln_dyninc      = .false. !  Logical switch for applying velocity increments
+    ln_sshinc      = .false. !  Logical switch for applying SSH increments
+    ln_asmdin      = .false. !  Logical switch for Direct Initialization (DI)
+    ln_asmiau      = .false. !  Logical switch for Incremental Analysis Updating (IAU)
+    ln_seaiceinc   = .false. !  Logical switch for applying sea ice increments
+    ln_phytobal    = .false. !  Logical switch for phytoplankton multivariate balancing
+    ln_slchltotinc = .false. !  Logical switch for applying slchltot increments
+    ln_slchldiainc = .false. !  Logical switch for applying slchldia increments
+    ln_slchlnoninc = .false. !  Logical switch for applying slchlnon increments
+    ln_schltotinc  = .false. !  Logical switch for applying schltot increments
+    ln_slphytotinc = .false. !  Logical switch for applying slphytot increments
+    ln_slphydiainc = .false. !  Logical switch for applying slphydia increments
+    ln_slphynoninc = .false. !  Logical switch for applying slphynon increments
+    ln_sfco2inc    = .false. !  Logical switch for applying sfCO2 increments
+    ln_spco2inc    = .false. !  Logical switch for applying spCO2 increments
+    ln_plchltotinc = .false. !  Logical switch for applying plchltot increments
+    ln_pchltotinc  = .false. !  Logical switch for applying pchltot increments
+    ln_pno3inc     = .false. !  Logical switch for applying pno3 increments
+    ln_psi4inc     = .false. !  Logical switch for applying psi4 increments
+    ln_pdicinc     = .false. !  Logical switch for applying pdic increments
+    ln_palkinc     = .false. !  Logical switch for applying palk increments
+    ln_pphinc      = .false. !  Logical switch for applying pph increments
+    ln_po2inc      = .false. !  Logical switch for applying po2 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]
-    nitdin    = 0          !  Timestep of background for DI in [0,nitend-nit000-1]
-    nitiaustr = 1          !  Timestep of start of IAU interval in [0,nitend-nit000-1]
-    nitiaufin = 15         !  Timestep of end of IAU interval in [0,nitend-nit000-1]
-    niaufn    = 0          !  Type of IAU weighting function
-    ln_salfix = .false.    !  Logical switch for ensuring that the sa > salfixmin
-    salfixmin = -9999      !  Minimum salinity after applying the increments
-    nn_divdmp = 0          !  Number of iterations of divergence damping operator
+    nitbkg         = 0       !  Timestep of background in [0,nitend-nit000-1]
+    nitdin         = 0       !  Timestep of background for DI in [0,nitend-nit000-1]
+    nitiaustr      = 1       !  Timestep of start of IAU interval in [0,nitend-nit000-1]
+    nitiaufin      = 15      !  Timestep of end of IAU interval in [0,nitend-nit000-1]
+    niaufn         = 0       !  Type of IAU weighting function
+    ln_salfix      = .false. !  Logical switch for ensuring that the sa > salfixmin
+    salfixmin      = -9999   !  Minimum salinity after applying the increments
+    nn_divdmp      = 0       !  Number of iterations of divergence damping operator
+    mld_choice_bgc = 1       !  MLD criterion to use for biogeochemistry assimilation
+    rn_maxchlinc   = -999.0  !  maximum absolute non-log chlorophyll increment from ocean colour assimilation
+                             !  <= 0 implies no maximum applied (switch turned off)
+                             !   > 0 implies maximum absolute chl increment capped at this value
 /
 !-----------------------------------------------------------------------
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 10149)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 10149)
@@ -0,0 +1,2195 @@
+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_roam'            : MEDUSA extras for carbonate cycle
+   !! 'key_karaml'          : Kara mixed layer depth
+   !!---------------------------------------------------------------------------
+   !! asm_bgc_check_options : check bgc assimilation options
+   !! asm_bgc_init_incs     : initialise bgc increments
+   !! asm_bgc_init_bkg      : initialise bgc background
+   !! asm_bgc_bal_wri       : write out bgc balancing increments
+   !! asm_bgc_bkg_wri       : write out bgc background
+   !! phyto2d_asm_inc       : apply the ocean colour increments
+   !! phyto3d_asm_inc       : apply the 3D phytoplankton increments
+   !! pco2_asm_inc          : apply the pCO2/fCO2 increments
+   !! ph_asm_inc            : apply the pH increments
+   !! bgc3d_asm_inc         : apply the generic 3D BGC increments
+   !!---------------------------------------------------------------------------
+   USE par_kind, ONLY:      & ! kind parameters
+      & wp
+   USE par_oce, ONLY:       & ! domain array sizes
+      & jpi,                &
+      & jpj,                &
+      & jpk
+   USE iom                    ! i/o
+   USE oce, ONLY:           & ! active tracer variables
+      & tsn
+   USE zdfmxl, ONLY :       & ! mixed layer depth 
+#if defined key_karaml
+      & hmld_kara,          &
+      & ln_kara,            &
+#endif   
+      & hmld,               & 
+      & hmlp,               &
+      & hmlpt 
+   USE asmpar, ONLY:        & ! assimilation parameters
+      & c_asmbkg,           &
+      & c_asmbal,           &
+      & nitbkg_r,           &
+      & nitdin_r,           &
+      & nitiaustr_r,        &
+      & nitiaufin_r
+#if defined key_top
+   USE trc, ONLY:           & ! passive tracer variables
+      & trn,                &
+      & trb
+   USE par_trc, ONLY:       & ! passive tracer parameters
+      & jptra
+#endif
+#if defined key_medusa
+   USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
+      & asm_phyto2d_bal_medusa
+   USE asmpco2bal, ONLY:    & ! pCO2 balancing for MEDUSA
+      & asm_pco2_bal
+   USE sms_medusa             ! MEDUSA parameters
+   USE par_medusa             ! MEDUSA parameters
+   USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system
+      & mocsy_interface
+   USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion
+      & f2pCO2
+   USE sms_medusa, ONLY:    & ! MEDUSA diagnostic variables
+      & pgrow_avg,          &
+      & ploss_avg,          &
+      & phyt_avg,           &
+      & mld_max
+#elif defined key_hadocc
+   USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC
+      & asm_phyto2d_bal_hadocc
+   USE asmpco2bal, ONLY:    & ! pCO2 balancing for HadOCC
+      & asm_pco2_bal
+   USE par_hadocc             ! HadOCC parameters
+   USE had_bgc_const          ! HadOCC parameters
+   USE trc, ONLY:           & ! HadOCC diagnostic variables
+      & pgrow_avg,          &
+      & ploss_avg,          &
+      & phyt_avg,           &
+      & mld_max,            &
+      & HADOCC_CHL
+   USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio
+      & cchl_p
+#endif
+
+   IMPLICIT NONE
+   PRIVATE                   
+
+   PUBLIC  asm_bgc_check_options  ! called by asm_inc_init in asminc.F90
+   PUBLIC  asm_bgc_init_incs      ! called by asm_inc_init in asminc.F90
+   PRIVATE asm_bgc_read_incs_2d   ! called by asm_bgc_init_incs
+   PRIVATE asm_bgc_read_incs_3d   ! called by asm_bgc_init_incs
+   PUBLIC  asm_bgc_init_bkg       ! called by asm_inc_init in asminc.F90
+   PUBLIC  asm_bgc_bal_wri        ! called by nemo_gcm in nemogcm.F90
+   PUBLIC  asm_bgc_bkg_wri        ! called by asm_bkg_wri in asmbkg.F90
+   PUBLIC  phyto2d_asm_inc        ! called by bgc_asm_inc in asminc.F90
+   PUBLIC  phyto3d_asm_inc        ! called by bgc_asm_inc in asminc.F90
+   PUBLIC  pco2_asm_inc           ! called by bgc_asm_inc in asminc.F90
+   PUBLIC  ph_asm_inc             ! called by bgc_asm_inc in asminc.F90
+   PUBLIC  bgc3d_asm_inc          ! called by bgc_asm_inc in asminc.F90
+
+   LOGICAL, PUBLIC :: ln_balwri      = .FALSE. !: No output of balancing incs
+   LOGICAL, PUBLIC :: ln_phytobal    = .FALSE. !: No phytoplankton balancing
+   LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No surface total      log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom     log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_schltotinc  = .FALSE. !: No surface total      chlorophyll        increment
+   LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total      log10(phyto C)     increment
+   LOGICAL, PUBLIC :: ln_slphydiainc = .FALSE. !: No surface diatom     log10(phyto C)     increment
+   LOGICAL, PUBLIC :: ln_slphynoninc = .FALSE. !: No surface non-diatom log10(phyto C)     increment
+   LOGICAL, PUBLIC :: ln_spco2inc    = .FALSE. !: No surface pCO2                          increment
+   LOGICAL, PUBLIC :: ln_sfco2inc    = .FALSE. !: No surface fCO2                          increment
+   LOGICAL, PUBLIC :: ln_plchltotinc = .FALSE. !: No profile total      log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_pchltotinc  = .FALSE. !: No profile total      chlorophyll        increment
+   LOGICAL, PUBLIC :: ln_pno3inc     = .FALSE. !: No profile nitrate                       increment
+   LOGICAL, PUBLIC :: ln_psi4inc     = .FALSE. !: No profile silicate                      increment
+   LOGICAL, PUBLIC :: ln_pdicinc     = .FALSE. !: No profile dissolved inorganic carbon    increment
+   LOGICAL, PUBLIC :: ln_palkinc     = .FALSE. !: No profile alkalinity                    increment
+   LOGICAL, PUBLIC :: ln_pphinc      = .FALSE. !: No profile pH                            increment
+   LOGICAL, PUBLIC :: ln_po2inc      = .FALSE. !: No profile oxygen                        increment
+
+   INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc 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]
+
+   REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc
+                                             ! <= 0 no maximum applied (switch turned off)
+                                             !  > 0 absolute chl inc capped at this value
+
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchltot_bkginc ! slchltot inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchldia_bkginc ! slchldia inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: schltot_bkginc  ! schltot inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphytot_bkginc ! slphytot inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphydia_bkginc ! slphydia inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphynon_bkginc ! slphynon inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: sfco2_bkginc    ! sfco2 inc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: spco2_bkginc    ! spco2 inc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: plchltot_bkginc ! plchltot inc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pchltot_bkginc  ! pchltot inc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pno3_bkginc     ! pno3 inc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: psi4_bkginc     ! psi4 inc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pdic_bkginc     ! pdic inc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: palk_bkginc     ! palk inc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pph_bkginc      ! pph inc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: po2_bkginc      ! po2 inc
+   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto2d_balinc  ! Balancing incs from ocean colour
+   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto3d_balinc  ! Balancing incs from p(l)chltot
+   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc     ! Balancing incs from spco2/sfco2
+   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ph_balinc       ! Balancing incs from pph
+
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg   ! Background phyto growth
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg   ! Background phyto loss
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg    ! Background phyto conc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg     ! Background max MLD
+   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg      ! Background tracer state
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: chl_bkg         ! Background chl
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: cchl_p_bkg      ! Background c:chl
+
+#  include "domzgr_substitute.h90"
+
+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 ! defined key_top || ( ! defined key_hadocc && ! defined key_medusa )
+      CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', &
+         &           ' but no compatible biogeochemical model is available' )
+#endif
+
+#if defined key_hadocc
+      IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slphydiainc .OR. &
+         & ln_slphynoninc .OR. ln_psi4inc .OR. ln_pphinc .OR. ln_po2inc ) THEN
+         CALL ctl_stop( ' Cannot assimilate PFTs, Si4, pH or O2 into HadOCC' )
+      ENDIF
+#endif
+
+#if defined key_medusa
+      IF ( .NOT. ln_foam_medusa ) THEN
+         CALL ctl_stop( ' MEDUSA data assimilation options not turned on: set ln_foam_medusa' )
+      ENDIF
+#endif
+
+      IF ( ( ln_phytobal ).AND.                                       &
+         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
+         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. &
+         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
+         & ( .NOT. ln_slphynoninc ) ) THEN
+         CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
+            &           ' if not assimilating ocean colour,',                   &
+            &           ' so ln_phytobal will be set to .false.')
+         ln_phytobal = .FALSE.
+      ENDIF
+
+      IF ( ( ln_balwri ).AND.                                         &
+         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
+         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. &
+         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
+         & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
+         & ( .NOT. ln_pchltotinc  ).AND.( .NOT. ln_pphinc      ).AND. &
+         & ( .NOT. ln_spco2inc    ).AND.( .NOT. ln_sfco2inc    ) ) THEN
+         CALL ctl_warn( ' Balancing increments not being calculated', &
+            &           ' so ln_balwri will be set to .false.')
+         ln_balwri = .FALSE.
+      ENDIF
+
+      IF ( ln_spco2inc .AND. ln_sfco2inc ) THEN
+         CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
+      ENDIF
+
+      IF ( ln_slchltotinc .AND. ln_schltotinc ) THEN
+         CALL ctl_stop( ' Can only assimilate surface log10(chlorophyll) or chlorophyll, not both' )
+      ENDIF
+
+      IF ( ln_plchltotinc .AND. ln_pchltotinc ) THEN
+         CALL ctl_stop( ' Can only assimilate profile log10(chlorophyll) or chlorophyll, not both' )
+      ENDIF
+
+      IF ( ( ln_slchltotinc .OR. ln_schltotinc  ) .AND. &
+         & ( ln_slchldiainc .OR. ln_slchlnoninc ) ) THEN
+         CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' )
+      ENDIF
+
+      IF ( ln_phytobal .AND.                                      &
+         & ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. &
+         &   ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN
+         CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
+            &           ' unless assimilating all model PFTs,')
+      ENDIF
+
+      IF ( ln_slphytotinc .AND. ( ln_slphydiainc .OR. ln_slphynoninc ) ) THEN
+         CALL ctl_stop( ' Can only assimilate total or PFT surface phytoplankton carbon, not both' )
+      ENDIF
+
+      IF ( ln_phytobal .AND.                                      &
+         & ( ( ln_slphynoninc .AND.( .NOT. ln_slphydiainc ) ).OR. &
+         &   ( ln_slphydiainc .AND.( .NOT. ln_slphynoninc ) ) ) ) THEN
+         CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
+            &           ' unless assimilating all model PFTs,')
+      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
+      !!
+      !!------------------------------------------------------------------------
+
+      ! Allocate and read increments
+      
+      IF ( ln_slchltotinc ) THEN
+         ALLOCATE( slchltot_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc )
+      ENDIF
+      
+      IF ( ln_slchldiainc ) THEN
+         ALLOCATE( slchldia_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc )
+      ENDIF
+      
+      IF ( ln_slchlnoninc ) THEN
+         ALLOCATE( slchlnon_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc )
+      ENDIF
+      
+      IF ( ln_schltotinc ) THEN
+         ALLOCATE( schltot_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc )
+      ENDIF
+      
+      IF ( ln_slphytotinc ) THEN
+         ALLOCATE( slphytot_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc )
+      ENDIF
+      
+      IF ( ln_slphydiainc ) THEN
+         ALLOCATE( slphydia_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc )
+      ENDIF
+      
+      IF ( ln_slphynoninc ) THEN
+         ALLOCATE( slphynon_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinslphynon', slphynon_bkginc )
+      ENDIF
+
+      IF ( ln_sfco2inc ) THEN
+         ALLOCATE( sfco2_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinsfco2', sfco2_bkginc )
+      ENDIF
+
+      IF ( ln_spco2inc ) THEN
+         ALLOCATE( spco2_bkginc(jpi,jpj) )
+         CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc )
+      ENDIF
+      
+      IF ( ln_plchltotinc ) THEN
+         ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) )
+         CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc )
+      ENDIF
+      
+      IF ( ln_pchltotinc ) THEN
+         ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) )
+         CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc )
+      ENDIF
+      
+      IF ( ln_pno3inc ) THEN
+         ALLOCATE( pno3_bkginc(jpi,jpj,jpk) )
+         CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc )
+      ENDIF
+      
+      IF ( ln_psi4inc ) THEN
+         ALLOCATE( psi4_bkginc(jpi,jpj,jpk) )
+         CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc )
+      ENDIF
+      
+      IF ( ln_pdicinc ) THEN
+         ALLOCATE( pdic_bkginc(jpi,jpj,jpk) )
+         CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc )
+      ENDIF
+      
+      IF ( ln_palkinc ) THEN
+         ALLOCATE( palk_bkginc(jpi,jpj,jpk) )
+         CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc )
+      ENDIF
+      
+      IF ( ln_pphinc ) THEN
+         ALLOCATE( pph_bkginc(jpi,jpj,jpk) )
+         CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc )
+      ENDIF
+      
+      IF ( ln_po2inc ) THEN
+         ALLOCATE( po2_bkginc(jpi,jpj,jpk) )
+         CALL asm_bgc_read_incs_3d( knum, 'bckinpo2', po2_bkginc )
+      ENDIF
+
+      ! Allocate balancing increments
+      
+      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
+         & ln_slphynoninc ) THEN
+#if defined key_top
+         ALLOCATE( phyto2d_balinc(jpi,jpj,jpk,jptra) )
+         phyto2d_balinc(:,:,:,:) = 0.0
+#else
+         CALL ctl_stop( ' key_top must be set for balancing increments' )
+#endif
+      ENDIF
+      
+      IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+#if defined key_top
+         ALLOCATE( phyto3d_balinc(jpi,jpj,jpk,jptra) )
+         phyto3d_balinc(:,:,:,:) = 0.0
+#else
+         CALL ctl_stop( ' key_top must be set for balancing increments' )
+#endif
+      ENDIF
+
+      IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
+#if defined key_top
+         ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
+         pco2_balinc(:,:,:,:) = 0.0
+#else
+         CALL ctl_stop( ' key_top must be set for balancing increments' )
+#endif
+      ENDIF
+
+      IF ( ln_pphinc ) THEN
+#if defined key_top
+         ALLOCATE( ph_balinc(jpi,jpj,jpk,jptra) )
+         ph_balinc(:,:,:,:) = 0.0
+#else
+         CALL ctl_stop( ' key_top must be set for balancing increments' )
+#endif
+      ENDIF
+
+   END SUBROUTINE asm_bgc_init_incs
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE asm_bgc_read_incs_2d( knum, cd_bgcname, p_incs )
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_init_incs  ***
+      !!
+      !! ** Purpose :   read 2d bgc increments
+      !!
+      !! ** Method  :   read increments from file
+      !!
+      !! ** Action  :   read increments from file
+      !!
+      !! References :   asm_inc_init
+      !!------------------------------------------------------------------------
+      !!
+      INTEGER,                      INTENT(in   ) :: knum       ! i/o unit
+      CHARACTER(LEN=13),            INTENT(in   ) :: cd_bgcname ! variable
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: p_incs     ! increments
+      !!
+      !!------------------------------------------------------------------------
+
+      ! Initialise
+      p_incs(:,:) = 0.0
+      
+      ! read from file
+      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 )
+      
+      ! Apply the masks
+      p_incs(:,:) = p_incs(:,:) * tmask(:,:,1)
+      
+      ! Set missing increments to 0.0 rather than 1e+20
+      ! to allow for differences in masks
+      WHERE( ABS( p_incs(:,:) ) > 1.0e+10 ) p_incs(:,:) = 0.0
+
+   END SUBROUTINE asm_bgc_read_incs_2d
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE asm_bgc_read_incs_3d( knum, cd_bgcname, p_incs )
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_init_incs  ***
+      !!
+      !! ** Purpose :   read 3d bgc increments
+      !!
+      !! ** Method  :   read increments from file
+      !!
+      !! ** Action  :   read increments from file
+      !!
+      !! References :   asm_inc_init
+      !!------------------------------------------------------------------------
+      !!
+      INTEGER,                          INTENT(in   ) :: knum       ! i/o unit
+      CHARACTER(LEN=13),                INTENT(in   ) :: cd_bgcname ! variable
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: p_incs     ! increments
+      !!
+      !!------------------------------------------------------------------------
+
+      ! Initialise
+      p_incs(:,:,:) = 0.0
+      
+      ! read from file
+      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 )
+      
+      ! Apply the masks
+      p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:)
+      
+      ! Set missing increments to 0.0 rather than 1e+20
+      ! to allow for differences in masks
+      WHERE( ABS( p_incs(:,:,:) ) > 1.0e+10 ) p_incs(:,:,:) = 0.0
+
+   END SUBROUTINE asm_bgc_read_incs_3d
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   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_top && ( defined key_hadocc || defined key_medusa )
+      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
+         & ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) 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,jpk)    )
+         ALLOCATE( cchl_p_bkg(jpi,jpj,jpk) )
+         chl_bkg(:,:,:)    = 0.0
+         cchl_p_bkg(:,:,:) = 0.0
+#endif
+         
+         !--------------------------------------------------------------------
+         ! Read background variables for phytoplankton 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(:,:,:)
+         cchl_p_bkg(:,:,:) = cchl_p_bkg(:,:,:) * tmask(:,:,:)
+#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) )
+         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) )
+#endif
+         
+         IF ( ln_phytobal ) 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_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 .OR. ln_pphinc ) 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 .OR. ln_pphinc ) 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
+#else
+      CALL ctl_stop( ' asm_bgc_init_bkg: key_top and a compatible biogeochemical model required' )
+#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 .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+            & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
+            & ln_slphynoninc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) )
+            CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) )
+            IF ( ln_phytobal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) )
+            ENDIF
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) )
+            IF ( ln_phytobal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) )
+               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
+            ENDIF
+#endif
+         ENDIF
+
+         IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) )
+            CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
+#endif
+         ENDIF
+
+         IF ( ln_spco2inc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ELSE IF ( ln_sfco2inc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ENDIF
+
+         IF ( ln_pphinc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) )
+#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(:,:,:)     )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         )
+#elif defined key_medusa
+      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg        )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg        )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg         )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max          )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd'  , trn(:,:,:,jpchd) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn'  , trn(:,:,:,jpphn) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd'  , trn(:,:,:,jpphd) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds'  , trn(:,:,:,jppds) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi'  , trn(:,:,:,jpzmi) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme'  , trn(:,:,:,jpzme) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din'  , trn(:,:,:,jpdin) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil'  , trn(:,:,:,jpsil) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer'  , trn(:,:,:,jpfer) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det'  , trn(:,:,:,jpdet) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc'  , trn(:,:,:,jpdtc) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic'  , trn(:,:,:,jpdic) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk'  , trn(:,:,:,jpalk) )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy'  , trn(:,:,:,jpoxy) )
+#endif
+
+   END SUBROUTINE asm_bgc_bkg_wri
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog )
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_init_incs  ***
+      !!
+      !! ** Purpose :   convert log increments to non-log
+      !!
+      !! ** Method  :   need to account for model background,
+      !!                cannot simply do 10^log_inc. Need to:
+      !!                1) Add log_inc to log10(background) to get log10(analysis)
+      !!                2) Take 10^log10(analysis) to get analysis
+      !!                3) Subtract background from analysis to get non-log incs
+      !!
+      !! ** Action  :   return non-log increments
+      !!
+      !! References :   
+      !!------------------------------------------------------------------------
+      !!
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pbkg        ! Background
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pinc_log    ! Log incs
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs
+      !
+      INTEGER                                     :: ji, jj      ! Loop counters
+      !!
+      !!------------------------------------------------------------------------
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( pbkg(ji,jj) > 0.0 ) THEN
+               pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + &
+                  &                       pinc_log(ji,jj) )      &
+                  &                 - pbkg(ji,jj)
+            ELSE
+               pinc_nonlog(ji,jj) = 0.0
+            ENDIF
+         END DO
+      END DO
+
+   END SUBROUTINE asm_bgc_unlog_2d
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE phyto2d_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
+      REAL(wp), DIMENSION(jpi,jpj) :: zmld        ! Mixed layer depth
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg
+#if defined key_medusa
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg
+      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs
+      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
+#endif
+      !!------------------------------------------------------------------------
+      
+      IF ( kt <= nit000 ) THEN
+
+         ! Un-log any log increments for passing to balancing routines
+         ! Remember that two sets of non-log increments should not be
+         ! expected to be in the same ratio as their log equivalents
+         
+         ! Total chlorophyll
+         IF ( ln_slchltotinc ) THEN
+#if defined key_medusa
+            zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
+#elif defined key_hadocc
+            zbkg_chltot(:,:) = chl_bkg(:,:,1)
+#endif
+            CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot )
+         ELSE IF ( ln_schltotinc ) THEN
+            zinc_chltot(:,:) = schltot_bkginc(:,:)
+         ELSE
+            zinc_chltot(:,:) = 0.0
+         ENDIF
+
+#if defined key_medusa
+         ! Diatom chlorophyll
+         IF ( ln_slchldiainc ) THEN
+            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
+            CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia )
+         ELSE
+            zinc_chldia(:,:) = 0.0
+         ENDIF
+#endif
+
+#if defined key_medusa
+         ! Non-diatom chlorophyll
+         IF ( ln_slchlnoninc ) THEN
+            zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn)
+            CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon )
+         ELSE
+            zinc_chlnon(:,:) = 0.0
+         ENDIF
+#endif
+
+         ! Total phytoplankton carbon
+         IF ( ln_slphytotinc ) THEN
+#if defined key_medusa
+            zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
+#elif defined key_hadocc
+            zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
+#endif
+            CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot )
+         ELSE
+            zinc_phytot(:,:) = 0.0
+         ENDIF
+
+#if defined key_medusa
+         ! Diatom phytoplankton carbon
+         IF ( ln_slphydiainc ) THEN
+            zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd
+            CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia )
+         ELSE
+            zinc_phydia(:,:) = 0.0
+         ENDIF
+#endif
+
+#if defined key_medusa
+         ! Non-diatom phytoplankton carbon
+         IF ( ln_slphynoninc ) THEN
+            zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn
+            CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon )
+         ELSE
+            zinc_phynon(:,:) = 0.0
+         ENDIF
+#endif
+
+         ! Select mixed layer
+         IF ( ll_asmdin ) THEN
+#if defined key_top && ( defined key_hadocc || defined key_medusa )
+            CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
+               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
+            zmld(:,:) = mld_max_bkg(:,:)
+#else
+            CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' )
+#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 phyto2d assimilation,', &
+                     &           ' but ln_kara=.false.' )
+               ENDIF
+#else
+               CALL ctl_stop( ' Kara mixed layer requested for phyto2d 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 phyto2d 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
+
+         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
+
+#if defined key_medusa
+         CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), &
+            &                         zinc_chltot,                         &
+            &                         ln_slchldiainc,                      &
+            &                         zinc_chldia,                         &
+            &                         ln_slchlnoninc,                      &
+            &                         zinc_chlnon,                         &
+            &                         ln_slphytotinc,                      &
+            &                         zinc_phytot,                         &
+            &                         ln_slphydiainc,                      &
+            &                         zinc_phydia,                         &
+            &                         ln_slphynoninc,                      &
+            &                         zinc_phynon,                         &
+            &                         zincper,                             &
+            &                         rn_maxchlinc, ln_phytobal, zmld,     &
+            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
+            &                         phyt_avg_bkg, mld_max_bkg,           &
+            &                         tracer_bkg, phyto2d_balinc )
+#elif defined key_hadocc
+         CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), &
+            &                         zinc_chltot,                         &
+            &                         ln_slphytotinc,                      &
+            &                         zinc_phytot,                         &
+            &                         zincper,                             &
+            &                         rn_maxchlinc, ln_phytobal, zmld,     &
+            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
+            &                         phyt_avg_bkg, mld_max_bkg,           &
+            &                         cchl_p_bkg(:,:,1),                   &
+            &                         tracer_bkg, phyto2d_balinc )
+#else
+         CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
+            &           '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,*) 'phyto2d_asm_inc : phyto2d 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
+            WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
+            END WHERE
+#elif defined key_hadocc
+            WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
+               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
+                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
+            END WHERE
+#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
+
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+#if defined key_medusa
+            WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
+               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
+               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+            END WHERE
+#elif defined key_hadocc
+            WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
+               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1)
+               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+            END WHERE
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE phyto2d_asm_inc
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE phyto3d_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the profile chlorophyll assimilation increments.
+      !!
+      !! ** Method  : Calculate increments to state variables.
+      !!              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                          :: ji, jj, jk   ! Loop counters
+      INTEGER                          :: it           ! Index
+      REAL(wp)                         :: zincwgt      ! IAU weight for timestep
+      REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn
+      REAL(wp)                         :: zfrac_chd    ! Fraction of jpchd
+      REAL(wp)                         :: zrat_phn_chn ! jpphn:jpchn ratio
+      REAL(wp)                         :: zrat_phd_chd ! jpphd:jpchd ratio
+      REAL(wp)                         :: zrat_pds_chd ! jppds:jpchd ratio
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc      ! Chlorophyll increments
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl      ! Chlorophyll background
+      !!------------------------------------------------------------------------
+
+      IF ( kt <= nit000 ) THEN
+
+         IF ( ln_plchltotinc ) THEN
+            ! Convert log10(chlorophyll) increment back to a chlorophyll increment
+            ! In order to transform logchl incs to chl incs, need to account for model
+            ! background, cannot simply do 10^logchl_bkginc. Need to:
+            ! 1) Add logchl inc to log10(background) to get log10(analysis)
+            ! 2) Take 10^log10(analysis) to get analysis
+            ! 3) Subtract background from analysis to get chl incs
+            ! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+#if defined key_medusa
+            bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
+#elif defined key_hadocc
+            bkg_chl(:,:,:) = chl_bkg(:,:,:)
+#endif
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
+                        chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
+                        IF ( rn_maxchlinc > 0.0 ) THEN
+                           chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
+                        ENDIF
+                     ELSE
+                        chl_inc(ji,jj,jk) = 0.0
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+         ELSE IF ( ln_pchltotinc ) THEN
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     IF ( rn_maxchlinc > 0.0 ) THEN
+                        chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
+                     ELSE
+                        chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
+                     ENDIF
+                  END DO
+               END DO
+            END DO
+         ENDIF
+
+#if defined key_medusa
+         ! Loop over each grid point partioning the increments based on existing ratios
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
+                     zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
+                     zfrac_chd = 1.0 - zfrac_chn
+                     phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
+                     phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
+                     zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
+                     zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
+                     zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
+                     phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
+                     phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
+                     phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
+                  ENDIF
+               END DO
+            END DO
+         END DO
+#elif defined key_hadocc
+         phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
+#else
+         CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
+            &           '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,*) 'phyto3d_asm_inc : phyto3d 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
+            WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
+            END WHERE
+#elif defined key_hadocc
+            WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
+               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
+                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
+            END WHERE
+#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
+
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+#if defined key_medusa
+            WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
+               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
+               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+            END WHERE
+#elif defined key_hadocc
+            WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
+               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1)
+               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+            END WHERE
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE phyto3d_asm_inc
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE pco2_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( sfco2_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(:,:) = sfco2_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 phytoplankton balancing if required
+         IF ( ln_phytobal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_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 phytoplankton balancing if required
+         IF ( ln_phytobal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_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 ( ll_asmdin ) THEN
+#if defined key_hadocc || defined key_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 ( 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
+            WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
+            END WHERE
+#elif defined key_hadocc
+            WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
+               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
+                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
+            END WHERE
+#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
+
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+#if defined key_medusa
+            WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
+               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)
+            END WHERE
+#elif defined key_hadocc
+            WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
+               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)
+            END WHERE
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE pco2_asm_inc
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
+      &                   ll_trainc, pt_bkginc, ps_bkginc )
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE ph_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the pH assimilation increments.
+      !!
+      !! ** Method  : Calculate increments to DIC and alkalinity from pH
+      !!              Use a similar approach to the pCO2 scheme
+      !!              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
+      
+      REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
+      REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
+      REAL(wp)                         :: PH_OUT         ! pH from pH calculation
+      REAL(wp)                         :: ph_bkg         ! pH from pH calculation
+      REAL(wp)                         :: ph_dic         ! pH from pH calculation
+      REAL(wp)                         :: ph_alk         ! pH from pH calculation
+      REAL(wp)                         :: dph_ddic       ! Change in pH wrt DIC
+      REAL(wp)                         :: dph_dalk       ! Change in pH wrt alk
+      REAL(wp)                         :: weight         ! Increment weighting
+      REAL(wp)                         :: zincwgt        ! IAU weight for current time step
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp   ! DIC background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp   ! Alkalinity background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp   ! DIN background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp   ! Silicate background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp   ! Temperature background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp   ! Salinity background
+      REAL(wp), DIMENSION(19)          :: dummy_out      ! Dummy variables
+      INTEGER                          :: ji, jj, jk, jx ! Loop counters
+      INTEGER                          :: it             ! Index
+      !!------------------------------------------------------------------------
+
+#if ! defined key_medusa
+      CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
+#else
+
+      IF ( kt <= nit000 ) THEN
+
+         !----------------------------------------------------------------------
+         ! DIC and alkalinity balancing
+         !----------------------------------------------------------------------
+
+         ! Account for physics assimilation if required
+         IF ( ll_trainc ) THEN
+            tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
+            sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
+         ELSE
+            tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
+            sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
+         ENDIF
+
+         ! Account for phytoplankton balancing if required
+         IF ( ln_phytobal ) THEN
+            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
+            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
+            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
+            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
+         ELSE
+            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
+            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
+            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
+            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
+         ENDIF
+         
+         ! Account for pCO2 balancing if required
+         IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
+            dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
+            alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
+         ENDIF
+         
+         ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
+         ! This requires three calls to the MOCSY carbonate package
+         ! One to find pH at background DIC and alkalinity
+         ! One to find pH when DIC is increased by 10
+         ! One to find pH when alkalinity is increased by 10
+         ! Then calculate the gradients
+
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
+
+                     DO jx = 1, 3
+                        IF ( jx == 1 ) THEN
+                           DIC_IN = dic_bkg_temp(ji,jj,jk)
+                           ALK_IN = alk_bkg_temp(ji,jj,jk)
+                        ELSE IF ( jx == 2 ) THEN
+                           DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
+                           ALK_IN = alk_bkg_temp(ji,jj,jk)
+                        ELSE IF ( jx == 3 ) THEN
+                           DIC_IN = dic_bkg_temp(ji,jj,jk)
+                           ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
+                        ENDIF
+
+                        CALL mocsy_interface(tsn(ji,jj,jk,jp_tem),         & ! IN:  temperature
+                           &                 tsn(ji,jj,jk,jp_sal),         & !      salinity
+                           &                 ALK_IN,                       & !      alkalinity
+                           &                 DIC_IN,                       & !      DIC
+                           &                 sil_bkg_temp(ji,jj,jk),       & !      silicate
+                           &                 din_bkg_temp(ji,jj,jk)/16.0,  & !      phosphate (Redfield from DIN)
+                           &                 1.0,                          & !      atmospheric pressure (dummy)
+                           &                 fsdept(ji,jj,jk),             & !      depth
+                           &                 gphit(ji,jj),                 & !      latitude
+                           &                 1.0,                          & !      gas transfer (dummy)
+                           &                 1.0,                          & !      atmospheric xCO2 (dummy)
+                           &                 1,                            & !      number of points
+                           &                 PH_OUT,                       & ! OUT: pH
+                           &                 dummy_out(1),  dummy_out(2),  & !      stuff we don't care about
+                           &                 dummy_out(3),  dummy_out(4),  &
+                           &                 dummy_out(5),  dummy_out(6),  &
+                           &                 dummy_out(7),  dummy_out(8),  &
+                           &                 dummy_out(9),  dummy_out(10), &
+                           &                 dummy_out(11), dummy_out(12), &
+                           &                 dummy_out(13), dummy_out(14), &
+                           &                 dummy_out(15), dummy_out(16), &
+                           &                 dummy_out(17), dummy_out(18), &
+                           &                 dummy_out(19) )
+
+                        IF ( jx == 1 ) THEN
+                           ph_bkg = PH_OUT
+                        ELSE IF ( jx == 2 ) THEN
+                           ph_dic = PH_OUT
+                        ELSE IF ( jx == 3 ) THEN
+                           ph_alk = PH_OUT
+                        ENDIF
+                     END DO
+
+                     dph_ddic = (ph_dic - ph_bkg) / zsearch
+                     dph_dalk = (ph_alk - ph_bkg) / zsearch
+                     weight   = pph_bkginc(ji,jj,jk) / ( (dph_ddic*dph_ddic) + (dph_dalk*dph_dalk) )
+
+                     ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
+                     ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
+                     
+                  ENDIF
+                  
+               END DO
+            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,*) 
+               WRITE(numout,*) 'ph_asm_inc : pH 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
+            WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
+               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
+            END WHERE
+
+            ! 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
+
+            ! 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 pH assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
+               &   trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
+               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1)
+               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+            END WHERE
+ 
+            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+#endif      
+      !
+   END SUBROUTINE ph_asm_inc
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
+   SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dyn_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
+      !!
+      !! ** Method  : 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
+
+         !----------------------------------------------------------------------
+         ! Remove any other balancing increments
+         !----------------------------------------------------------------------
+
+         IF ( ln_pno3inc ) THEN
+#if defined key_hadocc || defined key_medusa
+#if defined key_hadocc
+            it = jp_had_nut
+#elif defined key_medusa
+            it = jpdin
+#endif
+            IF ( ln_phytobal ) THEN
+               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
+               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_pphinc ) THEN
+               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it)
+            ENDIF
+#else
+            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+         ENDIF
+
+         IF ( ln_psi4inc ) THEN
+#if defined key_medusa
+            it = jpsil
+            IF ( ln_phytobal ) THEN
+               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
+               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_pphinc ) THEN
+               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
+            ENDIF
+#else
+            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+         ENDIF
+
+         IF ( ln_pdicinc ) THEN
+#if defined key_hadocc || defined key_medusa
+#if defined key_hadocc
+            it = jp_had_dic
+#elif defined key_medusa
+            it = jpdic
+#endif
+            IF ( ln_phytobal ) THEN
+               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
+               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_pphinc ) THEN
+               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it)
+            ENDIF
+#else
+            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+         ENDIF
+
+         IF ( ln_palkinc ) THEN
+#if defined key_hadocc || defined key_medusa
+#if defined key_hadocc
+            it = jp_had_alk
+#elif defined key_medusa
+            it = jpalk
+#endif
+            IF ( ln_phytobal ) THEN
+               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
+               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_pphinc ) THEN
+               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it)
+            ENDIF
+#else
+            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+         ENDIF
+
+         IF ( ln_po2inc ) THEN
+#if defined key_medusa
+            it = jpoxy
+            IF ( ln_phytobal ) THEN
+               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
+               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
+            ENDIF
+            IF ( ln_pphinc ) THEN
+               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it)
+            ENDIF
+#else
+            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+         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,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
+                  &  kt,' with IAU weight = ', pwgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the 3D BGC 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
+            ! Don't apply increments if they'll take concentrations negative
+
+            IF ( ln_pno3inc ) THEN
+#if defined key_hadocc
+               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
+               END WHERE
+#elif defined key_medusa
+               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+
+            IF ( ln_psi4inc ) THEN
+#if defined key_medusa
+               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+
+            IF ( ln_pdicinc ) THEN
+#if defined key_hadocc
+               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
+               END WHERE
+#elif defined key_medusa
+               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+
+            IF ( ln_palkinc ) THEN
+#if defined key_hadocc
+               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
+               END WHERE
+#elif defined key_medusa
+               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+
+            IF ( ln_po2inc ) THEN
+#if defined key_medusa
+               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
+                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
+                  trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+           
+            IF ( kt == nitiaufin_r ) THEN
+               IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
+               IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
+               IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
+               IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
+               IF ( ln_po2inc  ) DEALLOCATE( po2_bkginc  )
+            ENDIF
+
+         ENDIF
+
+      ELSEIF ( ll_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+#if defined key_hadocc
+            CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+#elif defined key_medusa
+            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+#endif
+
+            IF ( ln_pno3inc ) THEN
+#if defined key_hadocc
+               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
+                  trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
+               END WHERE
+#elif defined key_medusa
+               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
+                  trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+
+            IF ( ln_psi4inc ) THEN
+#if defined key_medusa
+               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
+                  trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+
+            IF ( ln_pdicinc ) THEN
+#if defined key_hadocc
+               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
+                  trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
+               END WHERE
+#elif defined key_medusa
+               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
+                  trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+
+            IF ( ln_palkinc ) THEN
+#if defined key_hadocc
+               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
+                  trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
+               END WHERE
+#elif defined key_medusa
+               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
+                  trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+
+            IF ( ln_po2inc ) THEN
+#if defined key_medusa
+               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
+                  &   trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
+                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
+                  trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
+               END WHERE
+#else
+               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
+#endif
+            ENDIF
+ 
+            IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
+            IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
+            IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
+            IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
+            IF ( ln_po2inc  ) DEALLOCATE( po2_bkginc  )
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE bgc3d_asm_inc
+
+   !!===========================================================================
+
+END MODULE asmbgc
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 10149)
@@ -49,4 +49,7 @@
 #if defined key_lim3
    USE ice
+#endif
+#if defined key_top
+   USE asmbgc, ONLY: asm_bgc_bkg_wri
 #endif
    USE timing
@@ -125,4 +128,8 @@
             IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
             !
+#if defined key_top
+            CALL asm_bgc_bkg_wri( kt, inum )
+            !
+#endif
             CALL iom_close( inum )
          ENDIF
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 10149)
@@ -44,4 +44,5 @@
 #endif
    USE sbc_oce          ! Surface boundary condition variables.
+   USE asmbgc           ! Biogeochemistry assimilation
 
    IMPLICIT NONE
@@ -54,4 +55,5 @@
    PUBLIC   ssh_asm_inc    !: Apply the SSH increment
    PUBLIC   seaice_asm_inc !: Apply the seaice increment
+   PUBLIC   bgc_asm_inc    !: Apply the biogeochemistry increments
 
 #if defined key_asminc
@@ -67,4 +69,5 @@
    LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
    LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration 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
@@ -133,10 +136,17 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace
       !!
-      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
+      NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri,                           &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
+         &                 ln_phytobal, ln_slchltotinc, ln_slchldiainc,    &
+         &                 ln_slchlnoninc, ln_schltotinc, ln_slphytotinc,  &
+         &                 ln_slphydiainc, ln_slphynoninc, ln_spco2inc,    &
+         &                 ln_sfco2inc, ln_plchltotinc, ln_pchltotinc,     &
+         &                 ln_pno3inc, ln_psi4inc, ln_pdicinc, ln_palkinc, &
+         &                 ln_pphinc, ln_po2inc,                           &
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
          &                 ln_salfix, salfixmin, nn_divdmp,                &
-         &                 ln_seaiceinc, ln_temnofreeze
+         &                 ln_seaiceinc, ln_temnofreeze,                   &
+         &                 mld_choice_bgc, rn_maxchlinc
       !!----------------------------------------------------------------------
 
@@ -161,6 +171,7 @@
          WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
+         WRITE(numout,*) '   Namelist nam_asminc : set assimilation increment parameters'
          WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
+         WRITE(numout,*) '      Logical switch for writing out balancing increments      ln_balwri = ', ln_balwri
          WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
          WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
@@ -168,4 +179,22 @@
          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 phytoplankton balancing             ln_phytobal = ', ln_phytobal
+         WRITE(numout,*) '      Logical switch for applying slchltot increments     ln_slchltotinc = ', ln_slchltotinc
+         WRITE(numout,*) '      Logical switch for applying slchldia increments     ln_slchldiainc = ', ln_slchldiainc
+         WRITE(numout,*) '      Logical switch for applying slchlnon increments     ln_slchlnoninc = ', ln_slchlnoninc
+         WRITE(numout,*) '      Logical switch for applying schltot increments       ln_schltotinc = ', ln_schltotinc
+         WRITE(numout,*) '      Logical switch for applying slphytot increments     ln_slphytotinc = ', ln_slphytotinc
+         WRITE(numout,*) '      Logical switch for applying slphydia increments     ln_slphydiainc = ', ln_slphydiainc
+         WRITE(numout,*) '      Logical switch for applying slphynon increments     ln_slphynoninc = ', ln_slphynoninc
+         WRITE(numout,*) '      Logical switch for applying spco2 increments           ln_spco2inc = ', ln_spco2inc
+         WRITE(numout,*) '      Logical switch for applying sfco2 increments           ln_sfco2inc = ', ln_sfco2inc
+         WRITE(numout,*) '      Logical switch for applying plchltot increments     ln_plchltotinc = ', ln_plchltotinc
+         WRITE(numout,*) '      Logical switch for applying pchltot increments       ln_pchltotinc = ', ln_pchltotinc
+         WRITE(numout,*) '      Logical switch for applying pno3 increments             ln_pno3inc = ', ln_pno3inc
+         WRITE(numout,*) '      Logical switch for applying psi4 increments             ln_psi4inc = ', ln_psi4inc
+         WRITE(numout,*) '      Logical switch for applying pdic increments             ln_pdicinc = ', ln_pdicinc
+         WRITE(numout,*) '      Logical switch for applying palk increments             ln_palkinc = ', ln_palkinc
+         WRITE(numout,*) '      Logical switch for applying pph increments               ln_pphinc = ', ln_pphinc
+         WRITE(numout,*) '      Logical switch for applying po2 increments               ln_po2inc = ', ln_po2inc
          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
@@ -176,4 +205,6 @@
          WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
          WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
+         WRITE(numout,*) '      Choice of MLD for BGC assimilation                  mld_choice_bgc = ', mld_choice_bgc
+         WRITE(numout,*) '      Maximum absolute chlorophyll increment (<=0 = off)    rn_maxchlinc = ', rn_maxchlinc
       ENDIF
 
@@ -212,4 +243,13 @@
          WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date
       ENDIF
+      
+      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
+         & ln_slphynoninc .OR. ln_spco2inc    .OR. ln_sfco2inc    .OR. &
+         & ln_plchltotinc .OR. ln_pchltotinc  .OR. ln_pno3inc     .OR. &
+         & ln_psi4inc     .OR. ln_pdicinc     .OR. ln_palkinc     .OR. &
+         & ln_pphinc      .OR. ln_po2inc ) THEN
+         lk_bgcinc = .TRUE.
+      ENDIF
 
       IF ( nacc /= 0 ) &
@@ -223,6 +263,8 @@
 
       IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
-           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
-         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
+         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. &
+         &        ( lk_bgcinc ) )) &
+         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' ln_(bgc-variable)inc is set to .true.', &
          &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
          &                ' Inconsistent options')
@@ -233,6 +275,7 @@
 
       IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
-         &                     )  &
-         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
+         & .AND.( .NOT. lk_bgcinc ) )  &
+         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' ln_(bgc-variable)inc are set to .false. :', &
          &                ' The assimilation increments are not applied')
 
@@ -254,4 +297,6 @@
          &                ' Background time step for Direct Initialization is outside', &
          &                ' the cycle interval')
+      
+      IF ( lk_bgcinc ) CALL asm_bgc_check_options
 
       IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
@@ -350,5 +395,6 @@
       ssh_iau(:,:)    = 0.0
 #endif
-      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN
+      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
+         &  .OR.( lk_bgcinc ) ) THEN
 
          !--------------------------------------------------------------------
@@ -424,4 +470,8 @@
          ENDIF
 
+         IF ( lk_bgcinc ) THEN
+            CALL asm_bgc_init_incs( inum )
+         ENDIF
+
          CALL iom_close( inum )
  
@@ -534,4 +584,8 @@
          CALL iom_close( inum )
 
+      ENDIF
+         
+      IF ( lk_bgcinc ) THEN
+         CALL asm_bgc_init_bkg
       ENDIF
       !
@@ -1155,4 +1209,57 @@
 
    END SUBROUTINE seaice_asm_inc
+
+
+   SUBROUTINE bgc_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE bgc_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 )
+      
+      ! Ocean colour variables first
+      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
+         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
+         & ln_slphynoninc ) THEN
+         CALL phyto2d_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau )
+      ENDIF
+      
+      ! Surface pCO2/fCO2 next
+      IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
+         CALL pco2_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, &
+            &               ln_trainc, t_bkginc, s_bkginc )
+      ENDIF
+      
+      ! Profile pH next
+      IF ( ln_pphinc ) THEN
+         CALL ph_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, &
+            &             ln_trainc, t_bkginc, s_bkginc )
+      ENDIF
+      
+      ! Then chlorophyll profiles
+      IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+         CALL phyto3d_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau )
+      ENDIF
+      
+      ! Remaining bgc profile variables
+      IF ( ln_pno3inc .OR. ln_psi4inc .OR. ln_pdicinc .OR. &
+         & ln_palkinc .OR. ln_po2inc ) THEN
+         CALL bgc3d_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau )
+      ENDIF
+
+   END SUBROUTINE bgc_asm_inc
    
    !!======================================================================
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90	(revision 10149)
@@ -20,7 +20,10 @@
       & c_asmtrj = 'assim_trj',                  & !: Filename for storing the 
                                                    !: reference trajectory
-      & c_asminc = 'assim_background_increments'   !: Filename for storing the 
+      & c_asminc = 'assim_background_increments', & !: Filename for storing the 
                                                    !: increments to the background
                                                    !: state
+      & c_asmbal = 'assim.balincs'                 !: Filename for storing the 
+                                                   !: balancing increments calculated
+                                                   !: for biogeochemistry
 
    INTEGER, PUBLIC :: nitbkg_r      !: Background time step referenced to nit000
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90	(revision 10149)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90	(revision 10149)
@@ -0,0 +1,366 @@
+MODULE asmpco2bal
+   !Contains subroutine to convert a global pco2 increment to independent DIC and alkalinity increments
+
+   USE par_kind, ONLY: wp
+   USE par_oce, ONLY: jpi, jpj, jpk
+   USE dom_oce, ONLY: tmask
+#if defined key_hadocc
+   USE had_bgc_const, ONLY: RHO_WATER_SI
+#else
+   USE phycst, ONLY: RHO_WATER_SI => rau0
+#endif
+   
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC asm_pco2_bal
+   CONTAINS 
+
+   SUBROUTINE  asm_pco2_bal(pco2inc, tco2_tracer_bkg, alk_tracer_bkg, &
+      &                     tem_tracer_bkg, sal_tracer_bkg,           &
+      &                     tco2_tracer_inc, alk_tracer_inc)
+      !Converts a global pco2 increment to independent DIC and alkalinity increments
+      !Calculates the solution to the equation [DDIC, DAlk]=(DpCO2 / |Grad (pCO2)|)* Grad(pCO2); where Grad(pco2)=[dpCO2/dDIC, dpCO2/dalk]
+      !Gradient code pco2_grad is a modification of assimcd1.mf77 by Ian Totterdell.
+
+      !James While 06/09
+
+      IMPLICIT NONE
+
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) ::  pco2inc            !pCO2 increment
+     
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg,      &    !In: tCO2 (DIC) background
+                                                  alk_tracer_bkg      !In alkalinity background
+     
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg,      &    !In: temperature background
+                                                  sal_tracer_bkg      !In salinity background
+     
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT) :: tco2_tracer_inc,      &    !Out: tCO2 (DIC) increment
+                                                   alk_tracer_inc      !Out alkalinity
+						   
+      !local variables
+      REAL(wp), DIMENSION(jpi,jpj) ::  Dpco2_tco2, Dpco2_alk
+      REAL(wp) :: wgt_fct
+      INTEGER :: ji,jj
+      
+      
+      !calculate the gradients with respect to DIC and Alkalinity
+      CALL pco2_grad(tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, Dpco2_tco2, Dpco2_alk)
+
+      
+      !loop over the ocean grid calculating the increments to tCO2 (DIC) and alkalinity
+      DO jj=1, jpj
+        DO ji=1, jpi
+
+          wgt_fct=pco2inc(ji,jj)/( (Dpco2_tco2(ji,jj)*Dpco2_tco2(ji,jj)) + (Dpco2_alk(ji,jj)*Dpco2_alk(ji,jj)) )
+          tco2_tracer_inc(ji,jj)=wgt_fct * Dpco2_tco2(ji,jj) * tmask(ji,jj,1)
+          alk_tracer_inc(ji,jj)=wgt_fct * Dpco2_alk(ji,jj) * tmask(ji,jj,1)
+ 
+        ENDDO
+      ENDDO
+       						  
+                                    
+
+   END SUBROUTINE  asm_pco2_bal
+   
+  
+   SUBROUTINE pco2_grad( tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, DPDC, DPDA )
+!
+!-------M---------M---------M---------M---------M---------M---------M--
+!
+!  Ian J Totterdell
+!  2007 August 08
+!
+!  Code for calculating derivatives of pCO2 wrt DIC and Alk
+!  
+!  J While 
+!  June 2009
+!  Modiefied for use in NEMO
+
+!-------M---------M---------M---------M---------M---------M---------M--
+      IMPLICIT NONE
+!
+!
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg,      &    !In: tCO2 (DIC) background
+                                                  alk_tracer_bkg      !In alkalinity background
+     
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg,      &    !In: temperature background
+                                                  sal_tracer_bkg      !In salinity background
+!
+      REAL(wp), INTENT(OUT), DIMENSION(jpi,jpj) :: &
+        DPDC,   &
+        DPDA
+!
+!
+!  Local variables.
+!
+      INTEGER :: &
+        mxiter,  &
+        ji,jj,   &
+        iter
+!
+      REAL(wp) :: &
+        a1_co2, a2_co2, a3_co2,               &
+        b1_co2, b2_co2, b3_co2,               &
+        a1_k1, a2_k1, a3_k1,                  &
+        b1_k1, b2_k1, b3_k1, b4_k1, b5_k1,    &
+        a1_k2, a2_k2, a3_k2,                  &
+        b1_k2, b2_k2, b3_k2, b4_k2, b5_k2,    &
+        a1_kb, a2_kb, a3_kb, a4_kb, a5_kb,    &
+        b1_kb, b2_kb, b3_kb,                  &
+        c1_kb, c2_kb, c3_kb, d_kb,            &
+        a1_kw, a2_kw, a3_kw,                  &
+        b1_kw, b2_kw, b3_kw, c_kw
+!
+      REAL(wp) :: &
+        tk, logtk, ttk,             &
+        sal, sqsal,                 &
+        kone, ktwo, kboric, kwater, &
+        alkprime, qa,qb,qc,         &
+        fxa, xtrma, xtrmb,          &
+        dxda, dxdc,                 &
+        dgda, dgdc
+!
+      REAL(wp), DIMENSION(jpi,jpj) :: &
+      co2_sol,     &
+      rek,         &
+      rg1,         &
+      rg2,         &
+      rg3,         &
+      boron,       &
+      tco2,        &
+      alk_t,       &
+      xs,          &
+      x1,          &
+      x2,          &
+      x2p,         &
+      x3,          &
+      xx,          &
+      xf
+!
+!  Set parameter and constant values.
+!
+      PARAMETER ( mxiter = 10 )
+!
+!  CO2 solubility, constants for the calculation of.
+!   From Weiss (1974)
+!
+      DATA &
+      a1_co2 / -60.2409 /,  a2_co2 / 93.4517 /,   a3_co2 / 23.3585 /, &
+      b1_co2 / 0.023517 /,  b2_co2 / -0.023656 /, b3_co2 / 0.0047036 /
+!
+!  kone, constants for the calculation of.
+!   From Roy et al (1993)
+!
+      DATA a1_k1 /-2307.1266/,  a2_k1 / 2.83655 /,   a3_k1 /-1.5529413/ &
+        , b1_k1 / -4.0484 /,   b2_k1 /-0.20760841/, b3_k1 /0.08468345/ &
+        , b4_k1 /-0.00654208/, b5_k1 / -0.001005 /
+!
+!  ktwo, constants for the calculation of.
+!   From Roy et al (1993)
+!
+      DATA a1_k2 /-3351.6106/,  a2_k2 / -9.226508 /, a3_k2 /-0.2005743/ &
+        , b1_k2 / -23.9722 /,  b2_k2 /-0.106901773/,b3_k2 /0.1130822/ &
+        , b4_k2 /-0.00846934/, b5_k2 / -0.001005 /
+!
+!  kboric, constants for the calculation of.
+!   Using data from Dickson (1990);
+!
+      DATA a1_kb / -8966.90 /,  a2_kb / -2890.53 /,  a3_kb / -77.942 / &
+        , a4_kb / 1.728 /,     a5_kb /-0.0996/ &
+        , b1_kb / 148.0248 /,  b2_kb / 137.1942 /,  b3_kb / 1.62142 / &
+        , c1_kb / -24.4344 /,  c2_kb / -25.085 /,   c3_kb / -0.2474 / &
+        , d_kb / 0.053105 /
+!
+!  kwater, constants for the calculation of.
+!   From Millero (1995, p670), using composite data.
+!
+      DATA a1_kw / -13847.26 /, a2_kw / 148.9652 /,  a3_kw / -23.6521 / &
+        , b1_kw / 118.67 /,    b2_kw / -5.977 /,    b3_kw / 1.0495 / &
+        , c_kw / -0.01615 / 
+!
+!  Executable code begins here.
+!
+!  Calculate the carbon system constants.
+!
+
+      DO ji=1,jpi
+        DO jj=1,jpj
+!
+          tk = 273.15 + MAX( -2.0, MIN( 40.0, tem_tracer_bkg(ji,jj) ) )
+          logtk = LOG(tk)
+          ttk = MAX( 2.71, 0.01 * ( 273.15 + tem_tracer_bkg(ji,jj) ) )
+          sal = MAX( 0.0, MIN( 45.0, sal_tracer_bkg(ji,jj) ) )
+          sqsal = SQRT(sal)
+!
+!  Calculate the CO2 solubility from Weiss (1974)
+!   Units are mol/(kg.soln)/atm .
+!   Check value (Zeebe and Wolf-G; tc=25,sal=35): ln(co2_sol) = -3.5617
+!
+          co2_sol(ji,jj) = EXP( &
+           a1_co2 + a2_co2 / ttk + a3_co2 * LOG(ttk) + &
+           sal * ( b1_co2 + ttk * ( b2_co2 + b3_co2 * ttk ) ) )
+!
+!  Calculate the 1st constant of carbonic acid from Roy et al. (1993)
+!   Units are mol/(kg.soln) on the total hydrogen ion scale.
+!   K1 = (HCO3-).(H+)/(CO2(g))
+!   Check value (Dickson and Goyet; tc=25,sal=35): ln(K1) = -13.4847
+!
+          kone &
+           = 1.0e-7 *( 1.0 + b5_k1*sal ) * EXP( a1_k1/tk + a2_k1 + a3_k1*logtk &
+           + ( b1_k1/tk + b2_k1 )*sqsal + b3_k1*sal + b4_k1*sal*sqsal )
+!
+!  Calculate the 2nd constant of carbonic acid from Roy et al. (1993)
+!   Units are mol/(kg.soln) on the total hydrogen ion scale.
+!   K2 = (CO3--).(H+)/(HCO3-)
+!   Check value (Dickson and Goyet; tc=25,sal=35): ln(K2) = -20.5504
+!
+          ktwo &
+           = 1.0e-7 *( 1.0 + b5_k2*sal ) * EXP( a1_k2/tk + a2_k2 + a3_k2*logtk &
+           + ( b1_k2/tk + b2_k2 )*sqsal + b3_k2*sal + b4_k2*sal*sqsal )
+!
+!  Calculate the equilibrium constant of boric acid from Dickson (1990).
+!   Units are mol/(kg.soln) on the total hydrogen ion scale.
+!   Kb = (B(OH)4-).(H+)/(B(OH)3) (with (Bt) = (B(OH)3) + (B(OH)4-) ).
+!   Check value (Dickson and Goyet; tc=25,sal=35): ln(Kb) = -19.7964
+!
+          kboric &
+           = 1.0e-7 * EXP( ( a1_kb + a2_kb*sqsal + a3_kb*sal + a4_kb*sal*sqsal &
+           + a5_kb*sal*sal )/tk + ( b1_kb + b2_kb*sqsal + b3_kb*sal ) &
+           + ( c1_kb + c2_kb*sqsal + c3_kb*sal )*logtk + d_kb*sqsal*tk )
+!
+!  Calculate the equilibrium constant of water from Millero (1994).
+!   Units are ( mol/(kg.soln) )^2 on the total hydrogen ion scale.
+!   Kw = (OH-).(H+)
+!   Check value (Dickson and Goyet; tc=25,sal=35): ln(Kw) = -30.434
+!
+          kwater &
+           = 1.0e-14  *EXP( ( a1_kw/tk + a2_kw + a3_kw*logtk ) &
+           + ( b1_kw/tk + b2_kw + b3_kw*logtk )*sqsal + c_kw*sal )
+!
+!  Calculate the combinations (as in PPCO2).
+!
+          rek(ji,jj) = SQRT( kone / ktwo )
+          rg3(ji,jj) = SQRT( kone * ktwo )
+          rg1(ji,jj) = rg3(ji,jj) / kboric
+          rg2(ji,jj) = kwater / rg3(ji,jj)
+!
+!  Set Boron concentration according to the surface salinity
+!   (Boron accounts for about 10% of the total Alkalinity).
+!   Use relationship used by Peng et al, Tellus 39b 1987,
+!   originally from Culkin 1965:
+!
+!         boron(ji,jj) = 4.16e-4 * sal / 35.0   ! Correct value but not in UM
+          boron(ji,jj) = 4.106e-4 * sal / 35.0
+!
+!  tco2 and alk_t are got from the modelled tracers at the surface
+!   (converted to units of mol/kg).
+!
+          tco2(ji,jj)  = MAX( 0.0, tco2_tracer_bkg(ji,jj) * 0.001 / RHO_WATER_SI )
+          alk_t(ji,jj) = MAX( 0.0, alk_tracer_bkg(ji,jj)  * 0.001 / RHO_WATER_SI )
+!
+!  Set initial values for iterative method.
+!
+          xs(ji,jj) = 1.0
+          x2p(ji,jj) = 0.0
+          x3(ji,jj) = 0.0
+          xx(ji,jj) = 0.0
+          xf(ji,jj) = xs(ji,jj)
+!
+        ENDDO
+      ENDDO
+
+!
+!  Use the secant method of similar triangles to calculate the H+ ion
+!   concentration, [H] (though strictly xf refers to  rg3/[H] ).
+!   To aid vectorisation this has been written to iterate a set
+!   number of times (mxiter), large enough to bring every example to
+!   within the 1e-5 fractional change criterion.
+!
+!  The method follows that of Bacastow (1990) in separating the
+!   alkalinity into two parts: the major part is the carbonate alk
+!   which uses the [H] that we are trying to find, while the minor
+!   part (xw) is sum of the borate and water dissociation terms which
+!   use the previous approximation to [H].
+!   This leaves an easy-to-solve quadratic.
+!
+      DO iter=1,mxiter
+!
+
+        DO ji=1,jpi
+          DO jj=1,jpj
+!
+            IF((abs(xf(ji,jj)-xx(ji,jj))/xf(ji,jj))<1.e-5) CYCLE   !prevent iteration collapsing to 0.0/0.0
+
+
+            x1(ji,jj) = x2p(ji,jj)
+            x2(ji,jj) = x3(ji,jj)
+            x2p(ji,jj) = xx(ji,jj)
+            x3(ji,jj) = xf(ji,jj)
+
+            IF (iter<=2) THEN
+              xx(ji,jj) = xf(ji,jj)
+            ELSE
+              xx(ji,jj) = x2p(ji,jj) + ( x3(ji,jj) - x2p(ji,jj) )*( x1(ji,jj) - x2p(ji,jj) ) &
+              /( x1(ji,jj) - x2(ji,jj) - x2p(ji,jj) + x3(ji,jj) )
+            ENDIF
+
+
+!
+!  All the above has been setting up the similar triangles for the
+!   secant method and using it to find an improved interpolated
+!   estimate; it has been completely independent of the equations
+!   we are trying to solve.
+!
+!  The next bit, however, IS dependent on the equations.
+!   alkprime is the carbonate alkalinity, which is calculated
+!   from the total alkalinity by subtracting the other, smaller, terms,
+!   which use the interpolated previous approximation to xf (xx).
+!
+
+            alkprime = alk_t(ji,jj) - rg2(ji,jj)*xx(ji,jj) + rg3(ji,jj)/xx(ji,jj) &
+            - boron(ji,jj)/( 1.0 + rg1(ji,jj)/xx(ji,jj) )	    
+	    
+!
+!  Set up the quadratic terms.
+!
+
+            qa = 2.0*tco2(ji,jj) - alkprime
+            qb = -rek(ji,jj)*( alkprime - tco2(ji,jj) )
+            qc = -alkprime
+
+            xf(ji,jj) = ( -qb + SQRT( qb*qb - 4.0*qa*qc ) )*0.5/qa
+
+
+          ENDDO
+        ENDDO
+!     
+      ENDDO
+
+!
+      DO ji=1,jpi
+        DO jj=1,jpj
+!
+          xtrma = 1.0/( 1.0 + rek(ji,jj)*xf(ji,jj) + xf(ji,jj)*xf(ji,jj) )
+          xtrmb = rg2(ji,jj) + rg3(ji,jj)/(xf(ji,jj)*xf(ji,jj)) &
+          + rg1(ji,jj)*boron(ji,jj)/(( rg1(ji,jj) + xf(ji,jj) )*( rg1(ji,jj) + xf(ji,jj) ))
+!
+          dxda = 1.0/( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) &
+          + rek(ji,jj)*xf(ji,jj)*xf(ji,jj) )*xtrma*xtrma + xtrmb )
+          dxdc = -dxda*xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma
+!
+          dgdc = xtrma - dxdc*tco2(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma*xtrma
+          dgda = ( 1.0 - dxda*( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) )*xtrma &
+          + xtrmb ) )/(xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) ))
+!
+          fxa = 1.0e-3*RHO_WATER_SI*co2_sol(ji,jj)
+          DPDC(ji,jj) = dgdc/fxa
+          DPDA(ji,jj) = dgda/fxa
+!
+        ENDDO
+      ENDDO
+         
+      RETURN
+      END SUBROUTINE pco2_grad
+
+END MODULE asmpco2bal
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_hadocc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_hadocc.F90	(revision 10149)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_hadocc.F90	(revision 10149)
@@ -0,0 +1,291 @@
+MODULE asmphyto2dbal_hadocc
+   !!======================================================================
+   !!                       ***  MODULE asmphyto2dbal_hadocc  ***
+   !! Calculate increments to HadOCC based on surface phyto2d increments
+   !!
+   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
+   !! For licensing reasons this is kept in its own internal Met Office
+   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
+   !! repository, and must be merged in when building.
+   !!
+   !!======================================================================
+   !! History :  3.6  ! 2017-08  (D. Ford)     Adapted from bioanal.F90
+   !!----------------------------------------------------------------------
+#if defined key_asminc && defined key_hadocc
+   !!----------------------------------------------------------------------
+   !! 'key_asminc'          : assimilation increment interface
+   !! 'key_hadocc'          : HadOCC model
+   !!----------------------------------------------------------------------
+   !! asm_phyto2d_bal_hadocc : routine to calculate increments to HadOCC
+   !!----------------------------------------------------------------------
+   USE par_kind,      ONLY: wp             ! kind parameters
+   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
+   USE dom_oce,       ONLY: gdepw_n        ! domain information
+   USE iom                                 ! i/o
+   USE par_hadocc                          ! HadOCC parameters
+   USE had_bgc_stnd,  ONLY: kmt            ! HadOCC parameters
+   USE had_bgc_const                       ! HadOCC parameters
+   USE par_trc,       ONLY: jptra          ! Tracer parameters
+   USE bioanalysis                         ! Nitrogen balancing
+
+   IMPLICIT NONE
+   PRIVATE                   
+
+   PUBLIC asm_phyto2d_bal_hadocc
+
+   ! Default values for biological assimilation parameters
+   ! Should match Hemmings et al. (2008)
+   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
+   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
+   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
+   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
+   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
+   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
+   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
+   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
+   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
+   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
+   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
+   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
+   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
+   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
+   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
+   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
+   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
+   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
+   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
+   !
+   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
+   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
+   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
+   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
+   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
+   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
+   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
+   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
+
+CONTAINS
+
+   SUBROUTINE asm_phyto2d_bal_hadocc( ld_chltot,                      &
+      &                               pinc_chltot,                    &
+      &                               ld_phytot,                      &
+      &                               pinc_phytot,                    &
+      &                               pincper,                        &
+      &                               p_maxchlinc, ld_phytobal, pmld, &
+      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                               phyt_avg_bkg, mld_max_bkg,      &
+      &                               cchl_p_bkg,                     &
+      &                               tracer_bkg, phyto2d_balinc )
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_phyto2d_bal_hadocc  ***
+      !!
+      !! ** Purpose :   calculate increments to HadOCC from 2d phytoplankton increments
+      !!
+      !! ** Method  :   call nitrogen balancing scheme
+      !!
+      !! ** Action  :   populate phyto2d_balinc
+      !!
+      !! References :   Hemmings et al., 2008, J. Mar. Res.
+      !!                Ford et al., 2012, Ocean Sci.
+      !!---------------------------------------------------------------------------
+      !!
+      LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments
+      LOGICAL,  INTENT(in   )                               :: ld_phytot      ! Assim phytot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phytot    ! phytot increments
+      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
+      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
+      LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: cchl_p_bkg     ! Surface C:Chl
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
+      !!
+      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
+      INTEGER                                               :: jkmax          ! Loop index
+      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
+      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
+      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
+      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
+      !!---------------------------------------------------------------------------
+
+      IF ( ( .NOT. ld_chltot ) .AND. ( .NOT. ld_phytot ) ) THEN
+         CALL ctl_stop( ' Trying to do phyto2d balancing but nothing to assimilate' )
+      ENDIF
+      
+      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      IF ( p_maxchlinc > 0.0 ) THEN
+         IF ( ld_chltot ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+
+      IF ( ld_phytot ) THEN
+         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
+      ENDIF
+      
+      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
+
+         ! Set up model parameters to be passed into Hemmings balancing routine
+         modparm(1)  = grow_sat
+         modparm(2)  = psmax
+         modparm(3)  = par
+         modparm(4)  = alpha
+         modparm(5)  = resp_rate
+         modparm(6)  = pmort_rate
+         modparm(7)  = phyto_min
+         modparm(8)  = z_mort_1
+         modparm(9)  = z_mort_2
+         modparm(10) = c2n_p
+         modparm(11) = c2n_z
+         modparm(12) = c2n_d
+         modparm(13) = graze_threshold
+         modparm(14) = holling_coef
+         modparm(15) = graze_sat
+         modparm(16) = graze_max
+
+         ! Set up assimilation parameters to be passed into balancing routine
+         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
+         assimparm(2)  = balnutext
+         assimparm(3)  = balnutmin
+         assimparm(4)  = r
+         assimparm(5)  = beta_g
+         assimparm(6)  = beta_l
+         assimparm(7)  = beta_m
+         assimparm(8)  = a_g
+         assimparm(9)  = a_l
+         assimparm(10) = a_m
+         assimparm(11) = zfracb0
+         assimparm(12) = zfracb1
+         assimparm(13) = qrfmax
+         assimparm(14) = qafmax
+         assimparm(15) = zrfmax
+         assimparm(16) = zafmax
+         assimparm(17) = prfmax
+         assimparm(18) = incphymin
+         assimparm(19) = integnstep
+         assimparm(20) = pthreshold
+
+         ! Set up external tracer indices array bstate
+         i_tracer(1) = 1   ! nutrient
+         i_tracer(2) = 2   ! phytoplankton
+         i_tracer(3) = 3   ! zooplankton
+         i_tracer(4) = 4   ! detritus
+         i_tracer(5) = 5   ! DIC
+         i_tracer(6) = 6   ! Alkalinity
+
+         ! Set background state
+         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_had_nut)
+         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_had_phy)
+         bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jp_had_zoo)
+         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_had_pdn)
+         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_had_dic)
+         bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jp_had_alk)
+
+         ! Call nitrogen balancing routine
+         CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,               &
+            &               n2be_p, n2be_z, n2be_d, assimparm,                          &
+            &               INT(pincper), 1, kmt(:,:), tmask(:,:,:),                    &
+            &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p_bkg(:,:), &
+            &               nbal_active, phyt_avg_bkg(:,:),                             &
+            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),          &
+            &               subsurf_active, deepneg_active,                             &
+            &               deeppos_active, nutprof_active,                             &
+            &               bstate, outincs,                                            &
+            &               diag_active, diag,                                          &
+            &               diag_fulldepth_active, diag_fulldepth )
+
+         ! Save balancing increments
+         phyto2d_balinc(:,:,:,jp_had_nut) = outincs(:,:,:,i_tracer(1))
+         phyto2d_balinc(:,:,:,jp_had_phy) = outincs(:,:,:,i_tracer(2))
+         phyto2d_balinc(:,:,:,jp_had_zoo) = outincs(:,:,:,i_tracer(3))
+         phyto2d_balinc(:,:,:,jp_had_pdn) = outincs(:,:,:,i_tracer(4))
+         phyto2d_balinc(:,:,:,jp_had_dic) = outincs(:,:,:,i_tracer(5))
+         phyto2d_balinc(:,:,:,jp_had_alk) = outincs(:,:,:,i_tracer(6))
+      
+      ELSE   ! No nitrogen balancing
+      
+         ! Initialise phytoplankton increment to zero
+         phyto2d_balinc(:,:,:,jp_had_phy) = 0.0
+         
+         ! Convert surface chlorophyll increment to phytoplankton nitrogen
+         phyto2d_balinc(:,:,1,jp_had_phy) = ( cchl_p_bkg(:,:) / (mw_carbon * c2n_p) ) * pinc_chltot(:,:)
+         
+         ! Propagate through mixed layer
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               jkmax = jpk-1
+               DO jk = jpk-1, 1, -1
+                  IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                     & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                     pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                     jkmax = jk
+                  ENDIF
+               END DO
+               !
+               DO jk = 2, jkmax
+                  phyto2d_balinc(ji,jj,jk,jp_had_phy) = phyto2d_balinc(ji,jj,1,jp_had_phy)
+               END DO
+               !
+            END DO
+         END DO
+
+         ! Set other balancing increments to zero
+         phyto2d_balinc(:,:,:,jp_had_nut) = 0.0
+         phyto2d_balinc(:,:,:,jp_had_zoo) = 0.0
+         phyto2d_balinc(:,:,:,jp_had_pdn) = 0.0
+         phyto2d_balinc(:,:,:,jp_had_dic) = 0.0
+         phyto2d_balinc(:,:,:,jp_had_alk) = 0.0
+      
+      ENDIF
+
+   END SUBROUTINE asm_phyto2d_bal_hadocc
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option : Empty routine
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE asm_phyto2d_bal_hadocc( ld_chltot,                      &
+      &                               pinc_chltot,                    &
+      &                               ld_phytot,                      &
+      &                               pinc_phytot,                    &
+      &                               pincper,                        &
+      &                               p_maxchlinc, ld_phytobal, pmld, &
+      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                               phyt_avg_bkg, mld_max_bkg,      &
+      &                               cchl_p_bkg,                     &
+      &                               tracer_bkg, phyto2d_balinc )
+      LOGICAL :: ld_chltot
+      REAL    :: pinc_chltot(:,:)
+      LOGICAL :: ld_phytot
+      REAL    :: pinc_phytot(:,:)
+      REAL    :: pincper
+      REAL    :: p_maxchlinc
+      LOGICAL :: ld_phytobal
+      REAL    :: pmld(:,:)
+      REAL    :: pgrow_avg_bkg(:,:)
+      REAL    :: ploss_avg_bkg(:,:)
+      REAL    :: phyt_avg_bkg(:,:)
+      REAL    :: mld_max_bkg(:,:)
+      REAL    :: cchl_p_bkg(:,:)
+      REAL    :: tracer_bkg(:,:,:,:)
+      REAL    :: phyto2d_balinc(:,:,:,:)
+      WRITE(*,*) 'asm_phyto2d_bal_hadocc: You should not have seen this print! error?'
+   END SUBROUTINE asm_phyto2d_bal_hadocc
+#endif
+
+   !!======================================================================
+END MODULE asmphyto2dbal_hadocc
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90	(revision 10149)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90	(revision 10149)
@@ -0,0 +1,510 @@
+MODULE asmphyto2dbal_medusa
+   !!======================================================================
+   !!                       ***  MODULE asmphyto2dbal_medusa  ***
+   !! Calculate increments to MEDUSA based on surface phyto2d increments
+   !!
+   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
+   !! For licensing reasons this is kept in its own internal Met Office
+   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
+   !! repository, and must be merged in when building.
+   !!
+   !!======================================================================
+   !! History :  3.6  ! 2017-08 (D. Ford)  Adapted from asmphyto2dbal_hadocc
+   !!----------------------------------------------------------------------
+#if defined key_asminc && defined key_medusa
+   !!----------------------------------------------------------------------
+   !! 'key_asminc'          : assimilation increment interface
+   !! 'key_medusa'          : MEDUSA model
+   !!----------------------------------------------------------------------
+   !! asm_phyto2d_bal_medusa : routine to calculate increments to MEDUSA
+   !!----------------------------------------------------------------------
+   USE par_kind,      ONLY: wp             ! kind parameters
+   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
+   USE dom_oce,       ONLY: gdepw_n        ! domain information
+   USE zdftmx,        ONLY: ln_tmx_itf, &  ! Indonesian Throughflow
+      &                     mask_itf       ! tidal mixing mask
+   USE iom                                 ! i/o
+   USE sms_medusa                          ! MEDUSA parameters
+   USE par_medusa                          ! MEDUSA parameters
+   USE par_trc,       ONLY: jptra          ! Tracer parameters
+   USE bioanalysis                         ! Nitrogen balancing
+
+   IMPLICIT NONE
+   PRIVATE                   
+
+   PUBLIC asm_phyto2d_bal_medusa
+
+   ! Default values for biological assimilation parameters
+   ! Should match Hemmings et al. (2008)
+   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
+   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
+   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
+   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
+   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
+   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
+   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
+   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
+   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
+   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
+   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
+   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
+   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
+   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
+   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
+   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
+   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
+   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
+   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
+   !
+   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
+   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
+   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
+   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
+   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
+   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
+   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
+   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
+
+CONTAINS
+
+   SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      &
+      &                               pinc_chltot,                    &
+      &                               ld_chldia,                      &
+      &                               pinc_chldia,                    &
+      &                               ld_chlnon,                      &
+      &                               pinc_chlnon,                    &
+      &                               ld_phytot,                      &
+      &                               pinc_phytot,                    &
+      &                               ld_phydia,                      &
+      &                               pinc_phydia,                    &
+      &                               ld_phynon,                      &
+      &                               pinc_phynon,                    &
+      &                               pincper,                        &
+      &                               p_maxchlinc, ld_phytobal, pmld, &
+      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                               phyt_avg_bkg, mld_max_bkg,      &
+      &                               tracer_bkg, phyto2d_balinc )
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_phyto2d_bal_medusa  ***
+      !!
+      !! ** Purpose :   calculate increments to MEDUSA from 2d phytoplankton increments
+      !!
+      !! ** Method  :   average up MEDUSA to look like HadOCC
+      !!                call nitrogen balancing scheme
+      !!                separate back out to MEDUSA
+      !!
+      !! ** Action  :   populate phyto2d_balinc
+      !!
+      !! References :   Hemmings et al., 2008, J. Mar. Res.
+      !!                Ford et al., 2012, Ocean Sci.
+      !!---------------------------------------------------------------------------
+      !!
+      LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments
+      LOGICAL,  INTENT(in   )                               :: ld_chldia      ! Assim chldia y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldia    ! chldia increments
+      LOGICAL,  INTENT(in   )                               :: ld_chlnon      ! Assim chlnon y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlnon    ! chlnon increments
+      LOGICAL,  INTENT(in   )                               :: ld_phytot      ! Assim phytot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phytot    ! phytot increments
+      LOGICAL,  INTENT(in   )                               :: ld_phydia      ! Assim phydia y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phydia    ! phydia increments
+      LOGICAL,  INTENT(in   )                               :: ld_phynon      ! Assim phynon y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phynon    ! phynon increments
+      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
+      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
+      LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
+      !!
+      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
+      INTEGER                                               :: jkmax          ! Loop index
+      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
+      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
+      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
+      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
+      REAL(wp)                                              :: zfrac          ! Fraction
+      REAL(wp)                                              :: zfrac_chn      ! Fraction of jpchn
+      REAL(wp)                                              :: zfrac_chd      ! Fraction of jpchd
+      REAL(wp)                                              :: zfrac_phn      ! Fraction of jpphn
+      REAL(wp)                                              :: zfrac_phd      ! Fraction of jpphd
+      REAL(wp)                                              :: zfrac_zmi      ! Fraction of jpzmi
+      REAL(wp)                                              :: zfrac_zme      ! Fraction of jpzme
+      REAL(wp)                                              :: zrat_pds_phd   ! Ratio of jppds:jpphd
+      REAL(wp)                                              :: zrat_chd_phd   ! Ratio of jpchd:jpphd
+      REAL(wp)                                              :: zrat_chn_phn   ! Ratio of jpchn:jpphn
+      REAL(wp)                                              :: zrat_phn_chn   ! Ratio of jpphn:jpchn
+      REAL(wp)                                              :: zrat_phd_chd   ! Ratio of jpphd:jpchd
+      REAL(wp)                                              :: zrat_pds_chd   ! Ratio of jppds:jpchd
+      REAL(wp)                                              :: zrat_dtc_det   ! Ratio of jpdtc:jpdet
+      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
+      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
+      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
+      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
+      !!---------------------------------------------------------------------------
+
+      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      IF ( p_maxchlinc > 0.0 ) THEN
+         IF ( ld_chltot ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
+               END DO
+            END DO
+         ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pinc_chltot(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnon(ji,jj)
+                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
+                  IF ( pinc_chltot(ji,jj) .NE. ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) ) THEN
+                     zfrac = pinc_chltot(ji,jj) / ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) )
+                     pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac
+                     pinc_chlnon(ji,jj) = pinc_chlnon(ji,jj) * zfrac
+                  ENDIF
+               END DO
+            END DO
+         ELSE IF ( ld_chldia ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pinc_chldia(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chldia(ji,jj), p_maxchlinc ) )
+                  pinc_chltot(ji,jj) = pinc_chldia(ji,jj)
+               END DO
+            END DO
+         ELSE IF ( ld_chlnon ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pinc_chlnon(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chlnon(ji,jj), p_maxchlinc ) )
+                  pinc_chltot(ji,jj) = pinc_chlnon(ji,jj)
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+
+      IF ( ld_phytot .OR. ld_phydia .OR. ld_phynon ) THEN
+         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
+      ENDIF
+      
+      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
+
+         ! Set up model parameters to be passed into Hemmings balancing routine.
+         ! For now these are hardwired to the standard HadOCC parameter values
+         ! (except C:N ratios) as this is what the scheme was developed for.
+         ! Obviously, HadOCC and MEDUSA are rather different models, so this
+         ! isn't ideal, but there's not always direct analogues between the two
+         ! parameter sets, so it's the easiest way to get something running.
+         ! In the longer term, some serious MarMOT-based development is required.
+         modparm(1)  = 0.1                               ! grow_sat
+         modparm(2)  = 2.0                               ! psmax
+         modparm(3)  = 0.845                             ! par
+         modparm(4)  = 0.02                              ! alpha
+         modparm(5)  = 0.05                              ! resp_rate
+         modparm(6)  = 0.05                              ! pmort_rate
+         modparm(7)  = 0.01                              ! phyto_min
+         modparm(8)  = 0.05                              ! z_mort_1
+         modparm(9)  = 1.0                               ! z_mort_2
+         modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p
+         modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z
+         modparm(12) = xthetad                           ! c2n_d
+         modparm(13) = 0.01                              ! graze_threshold
+         modparm(14) = 2.0                               ! holling_coef
+         modparm(15) = 0.5                               ! graze_sat
+         modparm(16) = 2.0                               ! graze_max
+
+         ! Set up assimilation parameters to be passed into balancing routine
+         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
+         assimparm(2)  = balnutext
+         assimparm(3)  = balnutmin
+         assimparm(4)  = r
+         assimparm(5)  = beta_g
+         assimparm(6)  = beta_l
+         assimparm(7)  = beta_m
+         assimparm(8)  = a_g
+         assimparm(9)  = a_l
+         assimparm(10) = a_m
+         assimparm(11) = zfracb0
+         assimparm(12) = zfracb1
+         assimparm(13) = qrfmax
+         assimparm(14) = qafmax
+         assimparm(15) = zrfmax
+         assimparm(16) = zafmax
+         assimparm(17) = prfmax
+         assimparm(18) = incphymin
+         assimparm(19) = integnstep
+         assimparm(20) = pthreshold
+
+         ! Set up external tracer indices array bstate
+         i_tracer(1) = 1   ! nutrient
+         i_tracer(2) = 2   ! phytoplankton
+         i_tracer(3) = 3   ! zooplankton
+         i_tracer(4) = 4   ! detritus
+         i_tracer(5) = 5   ! DIC
+         i_tracer(6) = 6   ! Alkalinity
+
+         ! Set background state
+         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin)
+         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd)
+         bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme)
+         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet)
+         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic)
+         bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk)
+
+         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
+         ! and nitrogen to biomass equivalent for PZD
+         ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA
+         cchl_p(:,:) = 0.0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN
+                  cchl_p(ji,jj) = xmassc * ( ( tracer_bkg(ji,jj,1,jpphn) * xthetapn ) +      &
+                     &                       ( tracer_bkg(ji,jj,1,jpphd) * xthetapd )   ) /  &
+                     &            ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) )
+               ENDIF
+            END DO
+         END DO
+         n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) )
+         n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) )
+         n2be_d = 14.01 + ( xmassc * xthetad )
+
+         ! Call nitrogen balancing routine
+         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
+            &               n2be_p, n2be_z, n2be_d, assimparm,                      &
+            &               INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
+            &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p(:,:), &
+            &               nbal_active, phyt_avg_bkg(:,:),                         &
+            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      &
+            &               subsurf_active, deepneg_active,                         &
+            &               deeppos_active, nutprof_active,                         &
+            &               bstate, outincs,                                        &
+            &               diag_active, diag,                                      &
+            &               diag_fulldepth_active, diag_fulldepth )
+         
+         ! Loop over each grid point partioning the increments
+         phyto2d_balinc(:,:,:,:) = 0.0
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  ! Phytoplankton
+                  IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. &
+                     & ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN
+                     IF ( ld_chltot ) THEN
+                        ! Phytoplankton nitrogen split up based on existing ratios
+                        zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / &
+                           &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
+                        zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / &
+                           &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
+                     ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
+                        ! Phytoplankton nitrogen split up based on assimilation increments
+                        zfrac_phn = pinc_chlnon(ji,jj) / pinc_chltot(ji,jj)
+                        zfrac_phd = pinc_chldia(ji,jj) / pinc_chltot(ji,jj)
+                     ENDIF
+
+                     ! Phytoplankton silicate split up based on existing ratios
+                     zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd)
+                     
+                     ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen
+                     ! Not using pinc_chltot directly as it's only 2D
+                     ! This method should give same results at surface as splitting pinc_chltot would
+                     zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn)
+                     zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd)
+                     
+                     phyto2d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn
+                     phyto2d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd
+                     phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
+                     phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
+                     phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_chd_phd
+                  ENDIF
+
+                  ! Zooplankton nitrogen split up based on existing ratios
+                  IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN
+                     zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / &
+                        &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
+                     zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / &
+                        &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
+                     phyto2d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi
+                     phyto2d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme
+                  ENDIF
+
+                  ! Nitrogen nutrient straight from balancing scheme
+                  phyto2d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1))
+
+                  ! Nitrogen detritus straight from balancing scheme
+                  phyto2d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4))
+
+                  ! DIC straight from balancing scheme
+                  phyto2d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5))
+
+                  ! Alkalinity straight from balancing scheme
+                  phyto2d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6))
+
+                  ! Remove diatom silicate increment from nutrient silicate to conserve mass
+                  IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto2d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
+                     phyto2d_balinc(ji,jj,jk,jpsil) = phyto2d_balinc(ji,jj,jk,jppds) * (-1.0)
+                  ENDIF
+
+                  ! Carbon detritus based on existing ratios
+                  IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN
+                     zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet)
+                     phyto2d_balinc(ji,jj,jk,jpdtc) = phyto2d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
+                  ENDIF
+
+                  ! Do nothing with iron or oxygen for the time being
+                  phyto2d_balinc(ji,jj,jk,jpfer) = 0.0
+                  phyto2d_balinc(ji,jj,jk,jpoxy) = 0.0
+                  
+               END DO
+            END DO
+         END DO
+      
+      ELSE   ! No nitrogen balancing
+      
+         ! Initialise individual chlorophyll increments to zero
+         phyto2d_balinc(:,:,:,jpchn) = 0.0
+         phyto2d_balinc(:,:,:,jpchd) = 0.0
+         
+         ! Split up total surface chlorophyll increments
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( ( tracer_bkg(ji,jj,1,jpchn) > 0.0 ) .AND. &
+                  & ( tracer_bkg(ji,jj,1,jpchd) > 0.0 ) ) THEN
+                  IF ( ld_chltot ) THEN
+                     ! Chlorophyll split up based on existing ratios
+                     zfrac_chn = tracer_bkg(ji,jj,1,jpchn) / &
+                        &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) )
+                     zfrac_chd = tracer_bkg(ji,jj,1,jpchd) / &
+                        &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) )
+                     phyto2d_balinc(ji,jj,1,jpchn) = pinc_chltot(ji,jj) * zfrac_chn
+                     phyto2d_balinc(ji,jj,1,jpchd) = pinc_chltot(ji,jj) * zfrac_chd
+                  ENDIF
+                  IF( ld_chldia ) THEN
+                     phyto2d_balinc(ji,jj,1,jpchd) = pinc_chldia(ji,jj)
+                  ENDIF
+                  IF( ld_chlnon ) THEN
+                     phyto2d_balinc(ji,jj,1,jpchn) = pinc_chlnon(ji,jj)
+                  ENDIF
+                  
+                  ! Maintain stoichiometric ratios of nitrogen and silicate
+                  IF ( ld_chltot .OR. ld_chlnon ) THEN
+                     zrat_phn_chn = tracer_bkg(ji,jj,1,jpphn) / tracer_bkg(ji,jj,1,jpchn)
+                     phyto2d_balinc(ji,jj,1,jpphn) = phyto2d_balinc(ji,jj,1,jpchn) * zrat_phn_chn
+                  ENDIF
+                  IF ( ld_chltot .OR. ld_chldia ) THEN
+                     zrat_phd_chd = tracer_bkg(ji,jj,1,jpphd) / tracer_bkg(ji,jj,1,jpchd)
+                     phyto2d_balinc(ji,jj,1,jpphd) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_phd_chd
+                     zrat_pds_chd = tracer_bkg(ji,jj,1,jppds) / tracer_bkg(ji,jj,1,jpchd)
+                     phyto2d_balinc(ji,jj,1,jppds) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_pds_chd
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+         
+         ! Propagate through mixed layer
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               jkmax = jpk-1
+               DO jk = jpk-1, 1, -1
+                  IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                     & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                     pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                     jkmax = jk
+                  ENDIF
+               END DO
+               !
+               DO jk = 2, jkmax
+                  phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,1,jpchn)
+                  phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,1,jpchd)
+                  phyto2d_balinc(ji,jj,jk,jpphn) = phyto2d_balinc(ji,jj,1,jpphn)
+                  phyto2d_balinc(ji,jj,jk,jpphd) = phyto2d_balinc(ji,jj,1,jpphd)
+                  phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,1,jppds)
+               END DO
+               !
+            END DO
+         END DO
+
+         ! Set other balancing increments to zero
+         phyto2d_balinc(:,:,:,jpzmi) = 0.0
+         phyto2d_balinc(:,:,:,jpzme) = 0.0
+         phyto2d_balinc(:,:,:,jpdin) = 0.0
+         phyto2d_balinc(:,:,:,jpsil) = 0.0
+         phyto2d_balinc(:,:,:,jpfer) = 0.0
+         phyto2d_balinc(:,:,:,jpdet) = 0.0
+         phyto2d_balinc(:,:,:,jpdtc) = 0.0
+         phyto2d_balinc(:,:,:,jpdic) = 0.0
+         phyto2d_balinc(:,:,:,jpalk) = 0.0
+         phyto2d_balinc(:,:,:,jpoxy) = 0.0
+
+      ENDIF
+      
+      ! If performing extra tidal mixing in the Indonesian Throughflow,
+      ! increments have been found to make the carbon cycle unstable
+      ! Therefore, mask these out
+      IF ( ln_tmx_itf ) THEN
+         DO jn = 1, jptra
+            DO jk = 1, jpk
+               phyto2d_balinc(:,:,jk,jn) = phyto2d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
+            END DO
+         END DO
+      ENDIF
+
+   END SUBROUTINE asm_phyto2d_bal_medusa
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option : Empty routine
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      &
+      &                              pinc_chltot,                    &
+      &                              ld_chldia,                      &
+      &                              pinc_chldia,                    &
+      &                              ld_chlnon,                      &
+      &                              pinc_chlnon,                    &
+      &                              ld_phytot,                      &
+      &                              pinc_phytot,                    &
+      &                              ld_phydia,                      &
+      &                              pinc_phydia,                    &
+      &                              ld_phynon,                      &
+      &                              pinc_phynon,                    &
+      &                              pincper,                        &
+      &                              p_maxchlinc, ld_phytobal, pmld, &
+      &                              pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                              phyt_avg_bkg, mld_max_bkg,      &
+      &                              tracer_bkg, phyto2d_balinc )
+      LOGICAL :: ld_chltot
+      REAL    :: pinc_chltot(:,:)
+      LOGICAL :: ld_chldia
+      REAL    :: pinc_chldia(:,:)
+      LOGICAL :: ld_chlnon
+      REAL    :: pinc_chlnon(:,:)
+      LOGICAL :: ld_phytot
+      REAL    :: pinc_phytot(:,:)
+      LOGICAL :: ld_phydia
+      REAL    :: pinc_phydia(:,:)
+      LOGICAL :: ld_phynon
+      REAL    :: pinc_phynon(:,:)
+      REAL    :: pincper
+      REAL    :: p_maxchlinc
+      LOGICAL :: ld_phytobal
+      REAL    :: pmld(:,:)
+      REAL    :: pgrow_avg_bkg(:,:)
+      REAL    :: ploss_avg_bkg(:,:)
+      REAL    :: phyt_avg_bkg(:,:)
+      REAL    :: mld_max_bkg(:,:)
+      REAL    :: tracer_bkg(:,:,:,:)
+      REAL    :: phyto2d_balinc(:,:,:,:)
+      WRITE(*,*) 'asm_phyto2d_bal_medusa: You should not have seen this print! error?'
+   END SUBROUTINE asm_phyto2d_bal_medusa
+#endif
+
+   !!======================================================================
+END MODULE asmphyto2dbal_medusa
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 10149)
@@ -38,15 +38,15 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: tidal mixing flag
 
-   !                       !!* Namelist  namzdf_tmx : tidal mixing *
-   REAL(wp) ::  rn_htmx     ! vertical decay scale for turbulence (meters)
-   REAL(wp) ::  rn_n2min    ! threshold of the Brunt-Vaisala frequency (s-1)
-   REAL(wp) ::  rn_tfe      ! tidal dissipation efficiency (St Laurent et al. 2002)
-   REAL(wp) ::  rn_me       ! mixing efficiency (Osborn 1980)
-   LOGICAL  ::  ln_tmx_itf  ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
-   REAL(wp) ::  rn_tfe_itf  ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   en_tmx     ! energy available for tidal mixing (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mask_itf   ! mask to use over Indonesian area
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   az_tmx     ! coefficient used to evaluate the tidal induced Kz
+   !                              !!* Namelist  namzdf_tmx : tidal mixing *
+   REAL(wp)        ::  rn_htmx     ! vertical decay scale for turbulence (meters)
+   REAL(wp)        ::  rn_n2min    ! threshold of the Brunt-Vaisala frequency (s-1)
+   REAL(wp)        ::  rn_tfe      ! tidal dissipation efficiency (St Laurent et al. 2002)
+   REAL(wp)        ::  rn_me       ! mixing efficiency (Osborn 1980)
+   LOGICAL, PUBLIC ::  ln_tmx_itf  ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
+   REAL(wp)        ::  rn_tfe_itf  ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   en_tmx     ! energy available for tidal mixing (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mask_itf   ! mask to use over Indonesian area
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)       ::   az_tmx     ! coefficient used to evaluate the tidal induced Kz
 
    !! * Substitutions
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 10149)
@@ -61,4 +61,5 @@
    USE asminc          ! assimilation increments     
    USE asmbkg          ! writing out state trajectory
+   USE asmbgc          ! biogeochemical assimilation increments
    USE diaptr          ! poleward transports           (dia_ptr_init routine)
    USE diadct          ! sections transports           (dia_dct_init routine)
@@ -160,4 +161,5 @@
                 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics
                 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
+                IF( lk_bgcinc ) CALL bgc_asm_inc( nit000 - 1 )    ! BGC
              ENDIF
           ENDIF
@@ -191,4 +193,6 @@
 
       IF( lk_diaobs   )   CALL dia_obs_wri
+      !
+      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/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 10149)
@@ -253,4 +253,6 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( lk_asminc .AND. ln_asmiau .AND. lk_bgcinc ) &
+         &               CALL bgc_asm_inc( kstp )     ! biogeochemistry assimilation
                          CALL trc_stp( kstp )         ! time-stepping
 #endif
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 10149)
@@ -74,7 +74,9 @@
                                    jriver_alk, jriver_c,                  &
                                    jriver_n, jriver_si,                   &
+                                   ln_foam_medusa,                        &
                                    riv_alk, riv_c, riv_n, riv_si,         &
                                    zn_dms_chd, zn_dms_chn, zn_dms_din,    &
                                    zn_dms_mld, zn_dms_qsr,                &
+                                   f2_pco2w, f2_fco2w,                    &
                                    xnln, xnld 
       USE trc,               ONLY: med_diag
@@ -87,4 +89,5 @@
 #  else
       USE trcco2_medusa,     ONLY: trc_co2_medusa
+      USE mocsy_mainmod,     ONLY: p2fCO2
 #  endif
       USE trcdms_medusa,     ONLY: trc_dms_medusa
@@ -327,4 +330,8 @@
                      iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt
                endif
+               IF ( ln_foam_medusa ) THEN
+                  !! DAF (Aug 2017): calculate fCO2 for observation operator
+                  CALL p2fCO2( f_pco2w, ztmp, f_pp0, 0.0, 1, f_fco2w )
+               ENDIF
             ENDIF
          ENDDO
@@ -506,4 +513,9 @@
                               CO2flux_conv
                !! ENDIF
+               IF ( ln_foam_medusa ) THEN
+                  !! DAF (Aug 2017): Save pCO2 and fCO2 for observation operator
+                  f2_pco2w(ji,jj) = f_pco2w(ji,jj)
+                  f2_fco2w(ji,jj) = f_pco2w(ji,jj)
+               ENDIF
                IF ( lk_iomput ) THEN
                   IF( med_diag%ATM_PCO2%dgsave ) THEN
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 10149)
@@ -46,4 +46,6 @@
                                    f3_co3, f3_h2co3, f3_hco3,           &
                                    f3_omarg, f3_omcal, f3_pH,           &
+                                   ln_foam_medusa, mld_max, pgrow_avg,  &
+                                   ploss_avg, phyt_avg,                 &
                                    za_sed_c, za_sed_ca, za_sed_fe,      &
                                    za_sed_n, za_sed_si,                 &
@@ -55,4 +57,5 @@
       USE trc,               ONLY: med_diag, nittrc000, trn 
       USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd
+      USE zdfmxl,            ONLY: hmld
  
       !! time (integer timestep)
@@ -113,4 +116,25 @@
 # endif      
 
+      IF ( ln_foam_medusa ) THEN
+         !!----------------------------------------------------------------------
+         !! Diagnostics required for ocean colour assimilation:
+         !! Mixed layer average phytoplankton growth, loss and concentration
+         !! Maximum mixed layer depth
+         !!----------------------------------------------------------------------
+         !!
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF ( hmld(ji,jj) .GT. 0.0 ) THEN
+                  pgrow_avg(ji,jj) = pgrow_avg(ji,jj) / hmld(ji,jj)
+                  ploss_avg(ji,jj) = ploss_avg(ji,jj) / hmld(ji,jj)
+                  phyt_avg(ji,jj)  = phyt_avg(ji,jj)  / hmld(ji,jj)
+                  IF ( hmld(ji,jj) .GT. mld_max(ji,jj) ) THEN
+                     mld_max(ji,jj) = hmld(ji,jj)
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+      ENDIF
+
 #  if defined key_debug_medusa
          !! AXY (12/07/17)
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 10149)
@@ -35,5 +35,5 @@
       USE bio_medusa_mod
       USE par_oce,           ONLY: jpi, jpj, jpk
-      USE sms_medusa,        ONLY: jdms
+      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max
       USE trc,               ONLY: ln_diatrc, med_diag, nittrc000 
       USE in_out_manager,    ONLY: lwp, numout
@@ -195,4 +195,11 @@
       fslowsinkc(:,:) = 0.0
 # endif      
+      !!
+      pgrow_avg(:,:) = 0.0
+      ploss_avg(:,:) = 0.0
+      phyt_avg(:,:)  = 0.0
+      IF( kt == nittrc000 ) THEN
+         mld_max(:,:) = 0.0
+      ENDIF
       !!
       !! allocate and initiate 2D diag
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 10149)
@@ -36,13 +36,22 @@
                                    fdpn, fdpn2, fdzme, fdzme2,             &
                                    fdzmi, fdzmi2, fsdiss, fsin,            &
+                                   fdep1, fprn, fprd,                      &
+                                   fgmepd, fgmepn, fgmipn,                 &
                                    zphd, zphn, zpds, zzme, zzmi
-      USE dom_oce,           ONLY: tmask
+      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask
+      USE par_kind,          ONLY: wp
       USE par_oce,           ONLY: jpim1, jpjm1
       USE phytoplankton_mod, ONLY: phytoplankton
       USE sms_medusa,        ONLY: jmpd, jmpn, jmzme, jmzmi,               &
+                                   ln_foam_medusa,                         &
+                                   pgrow_avg, ploss_avg, phyt_avg,         &
                                    xkphd, xkphn, xkzme, xkzmi,             &
                                    xmetapd, xmetapn, xmetazme, xmetazmi,   &
                                    xmpd, xmpn, xmzme, xmzmi, xsdiss
+      USE zdfmxl,            ONLY: hmld
       USE zooplankton_mod,   ONLY: zooplankton
+
+   !!* Substitution
+#  include "domzgr_substitute.h90"
 
       !! Level
@@ -50,4 +59,6 @@
 
       INTEGER :: ji, jj
+
+      REAL(wp) :: fq0
 
       !!-------------------------------------------------------------------
@@ -188,4 +199,39 @@
       ENDDO
 
+      IF ( ln_foam_medusa ) THEN
+         !! Mixed layer averages for ocean colour assimilation
+         !!
+         DO jj = 2,jpjm1
+            DO ji = 2,jpim1
+               IF (tmask(ji,jj,jk) == 1) THEN
+                  if (fdep1(ji,jj).le.hmld(ji,jj)) then
+                     !! this level is entirely in the mixed layer
+                     fq0 = 1.0
+                  elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then
+                     !! this level is entirely below the mixed layer
+                     fq0 = 0.0
+                  else
+                     !! this level straddles the mixed layer
+                     fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk)
+                  endif
+                  !!
+                  pgrow_avg(ji,jj) = pgrow_avg(ji,jj) +                   &
+                                     (((fprn(ji,jj) * zphn(ji,jj)) +      &
+                                       (fprd(ji,jj) * zphd(ji,jj))  ) *   &
+                                      fse3t(ji,jj,jk) * fq0)
+                  ploss_avg(ji,jj) = ploss_avg(ji,jj) +                   &
+                                     ((fgmepd(ji,jj) + fdpd(ji,jj) +      &
+                                       fdpd2(ji,jj)                +      &
+                                       fgmepn(ji,jj) + fdpn(ji,jj) +      &
+                                       fdpn2(ji,jj)  + fgmipn(ji,jj) ) *  &
+                                      fse3t(ji,jj,jk) * fq0)
+                  phyt_avg(ji,jj)  = phyt_avg(ji,jj)  +                   &
+                                     ((zphn(ji,jj) + zphd(ji,jj)) *       &
+                                      fse3t(ji,jj,jk) * fq0)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDIF
+
    END SUBROUTINE plankton
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 10149)
@@ -212,4 +212,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_ccd_arg  !: 2D aragonite CCD depth
 !!
+!! 2D fields of pCO2 and fCO2 for observation operator
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_pco2w    !: 2D pCO2
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_fco2w    !: 2D fCO2
+!!
 !! 2D fields of organic and inorganic material sedimented on the seafloor
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_sed_n    !: 2D organic nitrogen   (before)
@@ -289,4 +293,5 @@
    INTEGER  :: co2_rec
    REAL(wp) :: co2_yinit, co2_yend    !: First and Last year read in the xCO2.atm file
+   REAL(wp) :: xobs_xco2a             !: Observed atmospheric xCO2, from namelist
 #endif
 
@@ -349,4 +354,15 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: cmask       !: ???
 
+!!----------------------------------------------------------------------
+!! Parameters required for ocean colour assimilation
+!!----------------------------------------------------------------------
+!!
+   LOGICAL :: ln_foam_medusa                                 !: Flag to calculate and save diagnostics
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pgrow_avg  !: Mixed layer average phytoplankton growth
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ploss_avg  !: Mixed layer average phytoplankton loss
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: phyt_avg   !: Mixed layer average phytoplankton
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_max    !: Maximum mixed layer depth
+!!
+
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -361,5 +377,5 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_warn
-      INTEGER ::   ierr(8)        ! Local variables
+      INTEGER ::   ierr(9)        ! Local variables
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -371,4 +387,5 @@
       !* 2D and 3D fields of carbonate system parameters
       ALLOCATE( f2_ccd_cal(jpi,jpj)  , f2_ccd_arg(jpi,jpj)  ,       &
+                f2_pco2w(jpi,jpj)    , f2_fco2w(jpi,jpj)    ,       &
          &      f3_pH(jpi,jpj,jpk)   , f3_h2co3(jpi,jpj,jpk),       &
          &      f3_hco3(jpi,jpj,jpk) , f3_co3(jpi,jpj,jpk)  ,       &
@@ -419,4 +436,7 @@
          &      ffln(jpi,jpj,jpk)    , fflf(jpi,jpj,jpk)    ,       &
          &      ffls(jpi,jpj,jpk)    , cmask(jpi,jpj)       ,    STAT=ierr(8) ) 
+      !* Fields for ocean colour data assimilation
+      ALLOCATE( pgrow_avg(jpi,jpj)   , ploss_avg(jpi,jpj)   ,       &
+         &      phyt_avg(jpi,jpj)    , mld_max(jpi,jpj)     ,    STAT=ierr(9) )
 #endif
       !
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 10149)
@@ -105,5 +105,10 @@
       USE sbc_oce,                    ONLY: lk_oasis
       USE sms_medusa,                 ONLY: hist_pco2, co2_yinit, co2_yend, &
-                                            lk_pi_co2
+# if defined key_roam
+                                            xobs_xco2a,                     &
+# endif
+                                            pgrow_avg,                      &
+                                            ploss_avg, phyt_avg, mld_max,   &
+                                            lk_pi_co2, ln_foam_medusa
       USE trc,                        ONLY: ln_rsttr, nittrc000, trn
       USE bio_medusa_init_mod,        ONLY: bio_medusa_init
@@ -330,4 +335,7 @@
          !! f_xco2a(:,:) = 284.725       !! CMIP5 pre-industrial pCO2
          f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2 
+      ELSEIF ( xobs_xco2a > 0.0 ) THEN
+         IF(lwp) WRITE(numout,*) ' using observed atm pCO2 = ', xobs_xco2a
+         f_xco2a(:,:) = xobs_xco2a
       ELSE      
          !! xCO2 from file
@@ -418,6 +426,8 @@
       !! now use the NEMO calendar tool : nsec_month to be sure to call 
       !! at the beginning of a new month .
+      !! DAF: For FOAM we want to run daily
       IF ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
-           ( nsec_month .LE. INT(rdt) ) )  THEN
+           ( nsec_month .LE. INT(rdt) )          .OR.                        &
+           ( nsec_day   .LE. INT(rdt) .AND. ln_foam_medusa ) )  THEN
            IF ( lwp )  WRITE(numout,*)                                       &         
                               ' *** 3D carb chem call *** -- kt:', kt,       &
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90	(revision 10149)
@@ -330,5 +330,5 @@
       !!----------------------------------------------------------------------
       !! 
-      IF( ( .NOT.lk_oasis ) .AND. ( .NOT.lk_pi_co2 ) ) THEN
+      IF( ( .NOT.lk_oasis ) .AND. ( .NOT.lk_pi_co2 ) .AND. ( xobs_xco2a <= 0.0 ) ) THEN
          IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialisating atm CO2 record'
          CALL trc_ini_medusa_co2atm
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 10149)
@@ -87,7 +87,9 @@
       &  xsdiss,                                              &
       &  sedlam,sedlostpoc,jpkb,jdms,jdms_input,jdms_model,   &
-      &  scl_chl, chl_out, dmsmin, dmscut, dmsslp
+      &  scl_chl, chl_out, dmsmin, dmscut, dmsslp,            &
+      &  ln_foam_medusa
 #if defined key_roam
       NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        &
+      &    xobs_xco2a,                                        &
       &    xthetarem,xo2min 
 #endif
@@ -255,4 +257,6 @@
       dmscut      = 1.72 !! Anderson DMS default
       dmsslp      = 8.24 !! Anderson DMS default
+!!
+      ln_foam_medusa = .FALSE.
             
       !REWIND(numnatm)
@@ -415,4 +419,7 @@
 !!       chl_out     : select the chl field to send at the UM:
 !!                     1- Surf Chl ; 2- MLD Chl 
+!!
+!! FOAM - observation operator and data assimilation
+!!       ln_foam_medusa : calculate required diagnostics
 
       IF(lwp) THEN
@@ -975,4 +982,8 @@
             ENDIF
          ENDIF ! IF lk_oasis=true
+!! FOAM
+         WRITE(numout,*) '=== FOAM-related parameters'
+         WRITE(numout,*)     &
+         &   ' calculate diagnostics for data assimilation,            ln_foam_medusa = ', ln_foam_medusa
 !!
       ENDIF
@@ -1030,4 +1041,5 @@
       xthetarem = 0.
       xo2min    = 0.
+      xobs_xco2a = 0.
 
       !READ(numnatm,natroam)
@@ -1049,4 +1061,5 @@
 !!       xthetarem :  oxygen consumption by carbon remineralisation
 !!       xo2min    :  oxygen minimum concentration
+!!       xobs_xco2a : observed atmospheric xCO2 (not used if <= 0)
 
       IF(lwp) THEN
@@ -1066,4 +1079,6 @@
           WRITE(numout,*)     &
           &   ' oxygen minimum concentration                               xo2min      = ', xo2min
+          WRITE(numout,*)     &
+          &   ' observed atmospheric xCO2 (not used if <= 0)               xobs_xco2a  = ', xobs_xco2a
        ENDIF
 
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 10149)
@@ -36,5 +36,7 @@
    USE sbc_oce
    USE trcice          ! tracers in sea ice
+# if defined key_medusa
    USE sms_medusa      ! MEDUSA   initialisation
+# endif
    IMPLICIT NONE
    PRIVATE
@@ -107,5 +109,7 @@
       !                                                              ! total volume of the ocean 
       areatot = glob_sum( cvol(:,:,:) )
+# if defined key_medusa && defined key_roam
       carea(:,:) = e1e2t(:,:) * tmask(:,:,1) 
+# endif
 
       IF( lk_pisces  )       CALL trc_ini_pisces       ! PISCES  bio-model
Index: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 10020)
+++ branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 10149)
@@ -45,4 +45,5 @@
    USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable
    USE trcstat
+   USE obs_const, ONLY: obfillflt  ! Observation operator fill value
 
    IMPLICIT NONE
@@ -338,6 +339,34 @@
          IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...' 
       ENDIF
+      !
+      IF ( ln_foam_medusa ) THEN
+         !! 2D fields of pCO2 and fCO2 for observation operator on first timestep
+         IF( iom_varid( numrtr, 'PCO2W', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 present - reading in ...'
+            CALL iom_get( numrtr, jpdom_autoglo, 'PCO2W',  f2_pco2w(:,:)  )
+            CALL iom_get( numrtr, jpdom_autoglo, 'FCO2W',  f2_fco2w(:,:)  )
+         ELSE
+            IF(lwp) WRITE(numout,*) ' MEDUSA pCO2 absent - setting to fill ...'
+            f2_pco2w(:,:) = obfillflt * tmask(:,:,1)
+            f2_fco2w(:,:) = obfillflt * tmask(:,:,1)
+         ENDIF
+      ENDIF
 # endif
-
+      IF ( ln_foam_medusa ) THEN
+         !! Fields for ocean colour assimilation on first timestep
+         IF( iom_varid( numrtr, 'pgrow_avg', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg present - reading in ...'
+            CALL iom_get( numrtr, jpdom_autoglo, 'pgrow_avg',  pgrow_avg(:,:)  )
+            CALL iom_get( numrtr, jpdom_autoglo, 'ploss_avg',  ploss_avg(:,:)  )
+            CALL iom_get( numrtr, jpdom_autoglo, 'phyt_avg',   phyt_avg(:,:)   )
+            CALL iom_get( numrtr, jpdom_autoglo, 'mld_max',    mld_max(:,:)    )
+         ELSE
+            IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg absent - setting to zero ...'
+            pgrow_avg(:,:) = 0.0
+            ploss_avg(:,:) = 0.0
+            phyt_avg(:,:)  = 0.0
+            mld_max(:,:)   = 0.0
+         ENDIF
+      ENDIF
 
 #endif
@@ -510,5 +539,19 @@
       call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
       !!
+      IF ( ln_foam_medusa ) THEN
+         !! Fields for observation operator on first timestep
+         IF(lwp) WRITE(numout,*) ' MEDUSA OBS fields - writing out ...'
+         CALL iom_rstput( kt, nitrst, numrtw, 'PCO2W', f2_pco2w(:,:)  )
+         CALL iom_rstput( kt, nitrst, numrtw, 'FCO2W', f2_fco2w(:,:)  )
+      ENDIF
 # endif
+      IF ( ln_foam_medusa ) THEN
+         !! Fields for assimilation on first timestep
+         IF(lwp) WRITE(numout,*) ' MEDUSA ASM fields - writing out ...'
+         CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg', pgrow_avg(:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg', ploss_avg(:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg',  phyt_avg(:,:)  )
+         CALL iom_rstput( kt, nitrst, numrtw, 'mld_max',   mld_max(:,:)   )
+      ENDIF
 !!
 #endif
