Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9292)
@@ -1294,4 +1294,8 @@
     ln_asmiau = .false.    !  Logical switch for Incremental Analysis Updating (IAU)
     ln_seaiceinc = .false. !  Logical switch for applying sea ice increments
+    ln_logchlinc = .false. !  Logical switch for applying logchl increments
+    ln_logchlbal = .false. !  Logical switch for applying logchl multivariate balancing
+    ln_fco2inc = .false.   !  Logical switch for applying fCO2 increments
+    ln_pco2inc = .false.   !  Logical switch for applying pCO2 increments
     ln_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing
     nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1]
@@ -1303,4 +1307,8 @@
     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 logchl 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_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 9292)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 9292)
@@ -0,0 +1,134 @@
+MODULE asmbal
+   !!======================================================================
+   !!                       ***  MODULE asmbal  ***
+   !! Assimilation balancing interface: Write to file the balancing increments
+   !!                                   calculated for biogeochemistry
+   !!======================================================================
+   !!----------------------------------------------------------------------
+   !!   'key_asminc' : Switch on the assimilation increment interface
+   !!----------------------------------------------------------------------
+   !!   asm_bal_wri  : Write out the background state
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_kind, ONLY: &   ! Precision variables
+      & wp
+   USE iom                 ! I/O module
+   USE asminc              ! Main assimilation increments module
+   USE asmpar              ! Parameters for the assimilation interface
+#if defined key_medusa && defined key_foam_medusa
+   USE par_medusa
+#elif defined key_hadocc
+   USE par_hadocc
+#endif
+
+   IMPLICIT NONE
+
+   !! * Routine accessibility
+   PRIVATE
+   PUBLIC asm_bal_wri   !: Write out the background state
+
+CONTAINS
+
+   SUBROUTINE asm_bal_wri( kt )
+      !!-----------------------------------------------------------------------
+      !!
+      !!                  ***  ROUTINE asm_bal_wri ***
+      !!
+      !! ** Purpose : Write to file the balancing increments
+      !!              calculated for biogeochemistry
+      !!
+      !! ** Method  : Write to file the balancing increments
+      !!              calculated for biogeochemistry
+      !!
+      !! ** Action  :
+      !!                   
+      !! References :
+      !!
+      !! History    :
+      !!        ! 2014-08 (D. Ford)  Initial version, based on asm_bkg_wri
+      !!-----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( IN ) :: kt        ! Current time-step
+      !! * Local declarations
+      CHARACTER(LEN=50) :: cl_asmbal     ! Filename (with extension)
+      LOGICAL           :: llok          ! Check if file exists
+      INTEGER           :: inum          ! File unit number
+      REAL(wp)          :: zdate         ! Date
+      !!-----------------------------------------------------------------------
+     
+      ! Set things up
+      zdate = REAL( ndastp )
+      WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
+
+      ! Check if file exists
+      INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
+      IF( .NOT. llok ) THEN
+         IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
+            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
+
+         ! Define the output file       
+         CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
+
+         ! Write the information
+         CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate   )
+
+         IF ( ln_logchlinc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', logchl_balinc(:,:,:,jpchn) )
+            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', logchl_balinc(:,:,:,jpchd) )
+            IF ( ln_logchlbal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', logchl_balinc(:,:,:,jpphn) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', logchl_balinc(:,:,:,jpphd) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', logchl_balinc(:,:,:,jppds) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', logchl_balinc(:,:,:,jpzmi) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', logchl_balinc(:,:,:,jpzme) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', logchl_balinc(:,:,:,jpdin) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', logchl_balinc(:,:,:,jpsil) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', logchl_balinc(:,:,:,jpfer) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', logchl_balinc(:,:,:,jpdet) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', logchl_balinc(:,:,:,jpdtc) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', logchl_balinc(:,:,:,jpdic) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', logchl_balinc(:,:,:,jpalk) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', logchl_balinc(:,:,:,jpoxy) )
+            ENDIF
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', logchl_balinc(:,:,:,jp_had_phy) )
+            IF ( ln_logchlbal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', logchl_balinc(:,:,:,jp_had_nut) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', logchl_balinc(:,:,:,jp_had_zoo) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', logchl_balinc(:,:,:,jp_had_pdn) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', logchl_balinc(:,:,:,jp_had_dic) )
+               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', logchl_balinc(:,:,:,jp_had_alk) )
+            ENDIF
+#endif
+         ENDIF
+
+         IF ( ln_pco2inc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ELSE IF ( ln_fco2inc ) THEN
+#if defined key_medusa
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jpdic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jpalk) )
+#elif defined key_hadocc
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )
+#endif
+         ENDIF
+
+         CALL iom_close( inum )
+      ELSE
+         CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
+            &           ' Therefore not writing out balancing increments at this timestep', &
+            &           ' Something has probably gone wrong somewhere' )
+         IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
+            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
+      ENDIF
+                                 
+   END SUBROUTINE asm_bal_wri
+END MODULE asmbal
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 9292)
@@ -49,4 +49,21 @@
 #if defined key_lim3
    USE ice
+#endif
+#if defined key_hadocc
+   USE trc, ONLY: trn, &
+      &     pgrow_avg, &
+      &     ploss_avg, &
+      &     phyt_avg,  &
+      &     mld_max,   &
+      &     HADOCC_CHL
+   USE had_bgc_const, ONLY: cchl_p
+   USE par_hadocc
+#elif defined key_medusa && defined key_foam_medusa
+   USE trc, ONLY: trn
+   USE sms_medusa, ONLY: pgrow_avg, &
+      &                  ploss_avg, &
+      &                  phyt_avg,  &
+      &                  mld_max
+   USE par_medusa
 #endif
    IMPLICIT NONE
@@ -121,4 +138,38 @@
 !            CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               )
             CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
+#if defined key_hadocc
+            CALL iom_rstput( kt, nitbkg_r, inum, 'pgrow_avg'   , pgrow_avg             )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'ploss_avg'   , ploss_avg             )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'phyt_avg'    , phyt_avg              )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'mld_max'     , mld_max               )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_zoo'  , trn(:,:,:,jp_had_zoo) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_pdn'  , trn(:,:,:,jp_had_pdn) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_dic'  , trn(:,:,:,jp_had_dic) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_alk'  , trn(:,:,:,jp_had_alk) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_chl'  , HADOCC_CHL(:,:,1)     )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'hadocc_cchl' , cchl_p(:,:,1)         )
+#elif defined key_medusa && defined key_foam_medusa
+            CALL iom_rstput( kt, nitbkg_r, inum, 'pgrow_avg'   , pgrow_avg        )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'ploss_avg'   , ploss_avg        )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'phyt_avg'    , phyt_avg         )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'mld_max'     , mld_max          )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_chn'  , trn(:,:,:,jpchn) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_chd'  , trn(:,:,:,jpchd) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_phn'  , trn(:,:,:,jpphn) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_phd'  , trn(:,:,:,jpphd) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_pds'  , trn(:,:,:,jppds) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_zmi'  , trn(:,:,:,jpzmi) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_zme'  , trn(:,:,:,jpzme) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_din'  , trn(:,:,:,jpdin) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_sil'  , trn(:,:,:,jpsil) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_fer'  , trn(:,:,:,jpfer) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_det'  , trn(:,:,:,jpdet) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_dtc'  , trn(:,:,:,jpdtc) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_dic'  , trn(:,:,:,jpdic) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_alk'  , trn(:,:,:,jpalk) )
+            CALL iom_rstput( kt, nitbkg_r, inum, 'medusa_oxy'  , trn(:,:,:,jpoxy) )
+#endif
             !
             CALL iom_close( inum )
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 9292)
@@ -22,4 +22,6 @@
    !!   ssh_asm_inc    : Apply the SSH increment
    !!   seaice_asm_inc : Apply the seaice increment
+   !!   logchl_asm_inc : Apply the logchl increment
+   !!   pco2_asm_inc   : Apply the pCO2/fCO2 increment
    !!----------------------------------------------------------------------
    USE wrk_nemo         ! Memory Allocation
@@ -44,4 +46,34 @@
 #endif
    USE sbc_oce          ! Surface boundary condition variables.
+   USE zdfmxl, ONLY :  &  
+!   &  hmld_tref,       &   
+#if defined key_karaml
+   &  hmld_kara,       &
+   &  ln_kara,         &
+#endif   
+   &  hmld,            & 
+   &  hmlp,            &
+   &  hmlpt 
+#if defined key_top
+   USE trc, ONLY: &
+      & trn,      &
+      & trb
+   USE par_trc, ONLY: &
+      & jptra
+#endif
+#if defined key_medusa && defined key_foam_medusa
+   USE asmlogchlbal_medusa, ONLY: &
+      & asm_logchl_bal_medusa
+   USE asmpco2bal, ONLY: &
+      & asm_pco2_bal
+   USE par_medusa
+   USE mocsy_mainmod
+#elif defined key_hadocc
+   USE asmlogchlbal_hadocc, ONLY: &
+      & asm_logchl_bal_hadocc
+   USE asmpco2bal, ONLY: &
+      & asm_pco2_bal
+   USE par_hadocc
+#endif
 
    IMPLICIT NONE
@@ -54,4 +86,6 @@
    PUBLIC   ssh_asm_inc    !: Apply the SSH increment
    PUBLIC   seaice_asm_inc !: Apply the seaice increment
+   PUBLIC   logchl_asm_inc !: Apply the logchl increment
+   PUBLIC   pco2_asm_inc   !: Apply the pCO2/fCO2 increment
 
 #if defined key_asminc
@@ -61,4 +95,5 @@
 #endif
    LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields
+   LOGICAL, PUBLIC :: ln_balwri = .FALSE.      !: No output of the assimilation balancing increments
    LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment
    LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization
@@ -67,4 +102,8 @@
    LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
    LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment
+   LOGICAL, PUBLIC :: ln_logchlinc = .FALSE.   !: No log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_logchlbal = .FALSE.   !: Don't apply multivariate log10(chlorophyll) balancing
+   LOGICAL, PUBLIC :: ln_pco2inc = .FALSE.     !: No pCO2 increment
+   LOGICAL, PUBLIC :: ln_fco2inc = .FALSE.     !: No fCO2 increment
    LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
@@ -91,4 +130,32 @@
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
+
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pco2_bkginc   !: Increment to background pCO2
+#if defined key_top
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc    !: Increment to BGC variables from pCO2 assim
+#endif
+#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg  !: Background phyto growth for logchl balancing
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg  !: Background phyto loss for logchl balancing
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg   !: Background phyto for logchl balancing
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg    !: Background max MLD for logchl balancing
+   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg     !: Background tracer state variables
+#endif
+#if defined key_hadocc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: chl_bkg        !: Background surface chlorophyll
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: cchl_p_bkg     !: Background surface carbon:chlorophyll
+#endif
+   REAL(wp) :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation
+                                      !: <= 0 implies no maximum applied (switch turned off)
+                                      !:  > 0 implies maximum absolute chl increment capped at this value
+
+   INTEGER :: mld_choice_bgc    = 1   !: choice of mld criteria to use for biogeochemistry assimilation
+                                      !: 1) hmld      - Turbocline/mixing depth                           [W points]
+                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points]
+                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated]
+                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points]
+                                      !: 5) hmlpt     - Density criterion (0.01 kg/m^3 change from 10m)   [T points]
 
    !! * Substitutions
@@ -133,10 +200,13 @@
       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_logchlinc, ln_logchlbal,                     &
+         &                 ln_pco2inc, ln_fco2inc,                         &
          &                 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
       !!----------------------------------------------------------------------
 
@@ -163,4 +233,5 @@
          WRITE(numout,*) '   Namelist namasm : 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 +239,8 @@
          WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
          WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
+         WRITE(numout,*) '      Logical switch for applying logchl increments         ln_logchlinc = ', ln_logchlinc
+         WRITE(numout,*) '      Logical switch for logchl multivariate balancing      ln_logchlbal = ', ln_logchlbal
+         WRITE(numout,*) '      Logical switch for applying pCO2 increments             ln_pco2inc = ', ln_pco2inc
+         WRITE(numout,*) '      Logical switch for applying fCO2 increments             ln_fco2inc = ', ln_fco2inc
          WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
          WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
@@ -176,4 +251,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
 
@@ -223,6 +300,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. &
+         &        ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) )) &
+         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc is set to .true.', &
          &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
          &                ' Inconsistent options')
@@ -233,6 +312,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. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) )  &
+         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' ln_logchlinc, ln_pco2inc, and ln_fco2inc are set to .false. :', &
          &                ' The assimilation increments are not applied')
 
@@ -254,4 +334,20 @@
          &                ' Background time step for Direct Initialization is outside', &
          &                ' the cycle interval')
+
+      IF ( ( ln_balwri ).AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) ) THEN
+         CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', &
+            &           ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.')
+         ln_balwri = .FALSE.
+      ENDIF
+
+      IF ( ( ln_logchlbal ).AND.( .NOT. ln_logchlinc ) ) THEN
+         CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', &
+            &           ' Not assimilating logchl, so ln_logchlbal will be set to .false.')
+         ln_logchlbal = .FALSE.
+      ENDIF
+
+      IF ( ( ln_pco2inc ).AND.( ln_fco2inc ) ) THEN
+         CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
+      ENDIF
 
       IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
@@ -350,5 +446,22 @@
       ssh_iau(:,:)    = 0.0
 #endif
-      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN
+      IF ( ln_logchlinc ) THEN
+         ALLOCATE( logchl_bkginc(jpi,jpj) )
+         logchl_bkginc(:,:) = 0.0
+#if defined key_top
+         ALLOCATE( logchl_balinc(jpi,jpj,jpk,jptra) )
+         logchl_balinc(:,:,:,:) = 0.0
+#endif
+      ENDIF
+      IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
+         ALLOCATE( pco2_bkginc(jpi,jpj) )
+         pco2_bkginc(:,:) = 0.0
+#if defined key_top
+         ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
+         pco2_balinc(:,:,:,:) = 0.0
+#endif
+      ENDIF
+      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
+         &  .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN
 
          !--------------------------------------------------------------------
@@ -424,4 +537,26 @@
          ENDIF
 
+         IF ( ln_logchlinc ) THEN
+            CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl', logchl_bkginc(:,:), 1 )
+            ! Apply the masks
+            logchl_bkginc(:,:) = logchl_bkginc(:,:) * tmask(:,:,1)
+            ! Set missing increments to 0.0 rather than 1e+20
+            ! to allow for differences in masks
+            WHERE( ABS( logchl_bkginc(:,:) ) > 1.0e+10 ) logchl_bkginc(:,:) = 0.0
+         ENDIF
+
+         IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
+            IF ( ln_pco2inc ) THEN
+               CALL iom_get( inum, jpdom_autoglo, 'bckinpco2', pco2_bkginc(:,:), 1 )
+            ELSE IF ( ln_fco2inc ) THEN
+               CALL iom_get( inum, jpdom_autoglo, 'bckinfco2', pco2_bkginc(:,:), 1 )
+            ENDIF
+            ! Apply the masks
+            pco2_bkginc(:,:) = pco2_bkginc(:,:) * tmask(:,:,1)
+            ! Set missing increments to 0.0 rather than 1e+20
+            ! to allow for differences in masks
+            WHERE( ABS( pco2_bkginc(:,:) ) > 1.0e+10 ) pco2_bkginc(:,:) = 0.0
+         ENDIF
+
          CALL iom_close( inum )
  
@@ -535,4 +670,121 @@
 
       ENDIF
+
+#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
+      IF ( ln_logchlinc ) THEN
+
+         ALLOCATE( pgrow_avg_bkg(jpi,jpj)        )
+         ALLOCATE( ploss_avg_bkg(jpi,jpj)        )
+         ALLOCATE( phyt_avg_bkg(jpi,jpj)         )
+         ALLOCATE( mld_max_bkg(jpi,jpj)          )
+         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
+         pgrow_avg_bkg(:,:)  = 0.0
+         ploss_avg_bkg(:,:)  = 0.0
+         phyt_avg_bkg(:,:)   = 0.0
+         mld_max_bkg(:,:)    = 0.0
+         tracer_bkg(:,:,:,:) = 0.0
+
+#if defined key_hadocc
+         ALLOCATE( chl_bkg(jpi,jpj)    )
+         ALLOCATE( cchl_p_bkg(jpi,jpj) )
+         chl_bkg(:,:)    = 0.0
+         cchl_p_bkg(:,:) = 0.0
+#endif
+         
+         !--------------------------------------------------------------------
+         ! Read background variables for logchl assimilation
+         ! Some only required if performing balancing
+         !--------------------------------------------------------------------
+
+         CALL iom_open( c_asmbkg, inum )
+
+#if defined key_hadocc
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl',  chl_bkg    )
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
+         chl_bkg(:,:)    = chl_bkg(:,:)    * tmask(:,:,1)
+         cchl_p_bkg(:,:) = cchl_p_bkg(:,:) * tmask(:,:,1)
+#elif defined key_medusa
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
+#endif
+         
+         IF ( ln_logchlbal ) THEN
+
+            CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
+            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
+            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  )
+            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
+            pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
+            ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
+            phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1)
+            mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1)
+
+#if defined key_hadocc
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
+#elif defined key_medusa
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
+#endif
+         ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
+#if defined key_hadocc
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
+            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
+#elif defined key_medusa
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
+            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+#endif
+            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
+            mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
+         ENDIF
+
+         CALL iom_close( inum )
+         
+         DO jt = 1, jptra
+            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
+         END DO
+      
+      ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN
+
+         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
+         ALLOCATE( mld_max_bkg(jpi,jpj)          )
+         tracer_bkg(:,:,:,:) = 0.0
+         mld_max_bkg(:,:)    = 0.0
+
+         CALL iom_open( c_asmbkg, inum )
+         
+#if defined key_hadocc
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
+         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
+#elif defined key_medusa
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
+         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
+#endif
+         CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
+
+         CALL iom_close( inum )
+         
+         DO jt = 1, jptra
+            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
+         END DO
+         mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
+      
+      ENDIF
+#endif
       !
    END SUBROUTINE asm_inc_init
@@ -1155,4 +1407,382 @@
 
    END SUBROUTINE seaice_asm_inc
+
+   SUBROUTINE logchl_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE logchl_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the chlorophyll assimilation increments.
+      !!
+      !! ** Method  : Calculate increments to state variables using nitrogen
+      !!              balancing.
+      !!              Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(IN) :: kt   ! Current time step
+      !
+      INTEGER  :: jk              ! Loop counter
+      INTEGER  :: it              ! Index
+      REAL(wp) :: zincwgt         ! IAU weight for current time step
+      REAL(wp) :: zincper         ! IAU interval in seconds
+      !!----------------------------------------------------------------------
+
+      IF ( kt <= nit000 ) THEN
+
+         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
+
+#if defined key_medusa && defined key_foam_medusa
+         CALL asm_logchl_bal_medusa( logchl_bkginc, zincper, mld_choice_bgc, &
+            &                        rn_maxchlinc, ln_logchlbal, ln_asmdin,  &
+            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
+            &                        phyt_avg_bkg, mld_max_bkg,              &
+            &                        tracer_bkg, logchl_balinc )
+#elif defined key_hadocc
+         CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, &
+            &                        rn_maxchlinc, ln_logchlbal, ln_asmdin,  &
+            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
+            &                        phyt_avg_bkg, mld_max_bkg,              &
+            &                        chl_bkg, cchl_p_bkg,                    &
+            &                        tracer_bkg, logchl_balinc )
+#else
+         CALL ctl_stop( 'Attempting to assimilate logchl, ', &
+            &           'but not defined a biogeochemical model' )
+#endif
+
+      ENDIF
+
+      IF ( ln_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = wgtiau(it)      ! IAU weight for the current time step
+            ! note this is not a tendency so should not be divided by rdt
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
+                  &  kt,' with IAU weight = ', wgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the biogeochemical variables
+            ! Add directly to trn and trb, rather than to tra, because tra gets
+            ! reset to zero at the start of trc_stp, called after this routine
+#if defined key_medusa && defined key_foam_medusa
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+            END DO
+#elif defined key_hadocc
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
+                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
+                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+            END DO
+#endif
+
+            ! Do not deallocate arrays - needed by asm_bal_wri
+            ! which is called at end of model run
+         ENDIF
+
+      ELSEIF ( ln_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+#if defined key_medusa && defined key_foam_medusa
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+               &                         logchl_balinc(:,:,:,jp_msa0:jp_msa1)
+            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+#elif defined key_hadocc
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+               &                         logchl_balinc(:,:,:,jp_had0:jp_had1)
+            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE logchl_asm_inc
+
+
+   SUBROUTINE pco2_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE pco2_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
+      !!
+      !! ** Method  : Calculate increments to state variables using carbon
+      !!              balancing.
+      !!              Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(IN)                   :: kt           ! Current time step
+      !
+      INTEGER                               :: ji, jj, jk   ! Loop counters
+      INTEGER                               :: it           ! Index
+      INTEGER                               :: jkmax        ! Index of mixed layer depth
+      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
+      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
+      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
+      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
+      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
+      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
+      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
+      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
+      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
+
+      ! Coefficients for fCO2 to pCO2 conversion
+      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
+      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
+      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
+      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
+      REAL(wp)                              :: zcoef_fco2_5 = 2.0
+      REAL(wp)                              :: zcoef_fco2_6 = 57.7
+      REAL(wp)                              :: zcoef_fco2_7 = 0.118
+      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
+      !!----------------------------------------------------------------------
+
+      IF ( kt <= nit000 ) THEN
+
+         !--------------------------------------------------------------------
+         ! DIC and alkalinity balancing
+         !--------------------------------------------------------------------
+
+         IF ( ln_fco2inc ) THEN
+#if defined key_medusa && defined key_roam
+            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
+            patm(1) = 1.0
+            phyd(1) = 0.0
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  CALL f2pCO2( pco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
+               END DO
+            END DO
+#else
+            ! If assimilating fCO2, then convert to pCO2 using temperature
+            ! See flux_gas.F90 within HadOCC for details of calculation
+            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) /                                                             &
+               &                    EXP((zcoef_fco2_1                                                            + &
+               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
+               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
+               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
+               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
+               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
+#endif
+         ELSE
+            pco2_bkginc_temp(:,:) = pco2_bkginc(:,:)
+         ENDIF
+
+         ! Account for physics assimilation if required
+         IF ( ln_trainc ) THEN
+            tem_bkg_temp(:,:) = tsn(:,:,1,1) + t_bkginc(:,:,1)
+            sal_bkg_temp(:,:) = tsn(:,:,1,2) + s_bkginc(:,:,1)
+         ELSE
+            tem_bkg_temp(:,:) = tsn(:,:,1,1)
+            sal_bkg_temp(:,:) = tsn(:,:,1,2)
+         ENDIF
+
+#if defined key_medusa
+         ! Account for logchl balancing if required
+         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + logchl_balinc(:,:,1,jpdic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + logchl_balinc(:,:,1,jpalk)
+         ELSE
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
+         ENDIF
+
+         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
+            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
+            &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
+
+#elif defined key_hadocc
+         ! Account for logchl balancing if required
+         IF ( ln_logchlinc .AND. ln_logchlbal ) THEN
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + logchl_balinc(:,:,1,jp_had_dic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + logchl_balinc(:,:,1,jp_had_alk)
+         ELSE
+            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
+            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
+         ENDIF
+
+         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
+            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
+            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
+
+#else
+         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
+            &           'but not defined a biogeochemical model' )
+#endif
+
+         ! Select mixed layer
+         IF ( ln_asmdin ) THEN
+#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
+            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
+               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
+            zmld(:,:) = mld_max_bkg(:,:)
+#else
+            CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
+#endif
+         ELSE
+            SELECT CASE( mld_choice_bgc )
+            CASE ( 1 )                   ! Turbocline/mixing depth [W points]
+               zmld(:,:) = hmld(:,:)
+            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
+               zmld(:,:) = hmlp(:,:)
+            CASE ( 3 )                   ! Kara MLD [Interpolated]
+#if defined key_karaml
+               IF ( ln_kara ) THEN
+                  zmld(:,:) = hmld_kara(:,:)
+               ELSE
+                  CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
+                     &           ' but ln_kara=.false.' )
+               ENDIF
+#else
+               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
+                  &           ' but is not defined' )
+#endif
+            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
+               !zmld(:,:) = hmld_tref(:,:)
+               CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
+                  &           ' but is not available in this version' )
+            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
+               zmld(:,:) = hmlpt(:,:)
+            END SELECT
+         ENDIF
+
+         ! Propagate through mixed layer
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               jkmax = jpk-1
+               DO jk = jpk-1, 1, -1
+                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                     jkmax = jk
+                  ENDIF
+               END DO
+               !
+#if defined key_top
+               DO jk = 2, jkmax
+                  pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
+               END DO
+#endif
+               !
+            END DO
+         END DO
+
+      ENDIF
+
+      IF ( ln_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = wgtiau(it)      ! IAU weight for the current time step
+            ! note this is not a tendency so should not be divided by rdt
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               IF ( ln_pco2inc ) THEN
+                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
+                     &  kt,' with IAU weight = ', wgtiau(it)
+               ELSE IF ( ln_fco2inc ) THEN
+                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
+                     &  kt,' with IAU weight = ', wgtiau(it)
+               ENDIF
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the biogeochemical variables
+            ! Add directly to trn and trb, rather than to tra, because tra gets
+            ! reset to zero at the start of trc_stp, called after this routine
+#if defined key_medusa && defined key_foam_medusa
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+            END DO
+#elif defined key_hadocc
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
+                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
+                  &                          pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+            END DO
+#endif
+
+            ! Do not deallocate arrays - needed by asm_bal_wri
+            ! which is called at end of model run
+
+         ENDIF
+
+      ELSEIF ( ln_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+#if defined key_medusa && defined key_foam_medusa
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+               &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1)
+            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+#elif defined key_hadocc
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+               &                         pco2_balinc(:,:,:,jp_had0:jp_had1)
+            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE pco2_asm_inc
    
    !!======================================================================
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_hadocc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_hadocc.F90	(revision 9292)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_hadocc.F90	(revision 9292)
@@ -0,0 +1,317 @@
+MODULE asmlogchlbal_hadocc
+   !!======================================================================
+   !!                       ***  MODULE asmlogchlbal_hadocc  ***
+   !! Calculate increments to HadOCC based on surface logchl 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_logchl_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 zdfmxl                              ! mixed layer depth
+   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_logchl_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_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, &
+      &                              k_maxchlinc, ld_logchlbal, ld_asmdin,   &
+      &                              pgrow_avg_bkg, ploss_avg_bkg,           &
+      &                              phyt_avg_bkg, mld_max_bkg,              &
+      &                              chl_bkg, cchl_p_bkg,                    &
+      &                              tracer_bkg, logchl_balinc )
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_logchl_bal_hadocc  ***
+      !!
+      !! ** Purpose :   calculate increments to HadOCC from logchl increments
+      !!
+      !! ** Method  :   convert logchl increments to chl increments
+      !!                call nitrogen balancing scheme
+      !!
+      !! ** Action  :   populate logchl_balinc
+      !!
+      !! References :   Hemmings et al., 2008, J. Mar. Res.
+      !!                Ford et al., 2012, Ocean Sci.
+      !!---------------------------------------------------------------------------
+      !!
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: logchl_bkginc  ! logchl increments
+      REAL(wp), INTENT(in   )                               :: aincper        ! Assimilation period
+      INTEGER,  INTENT(in   )                               :: mld_choice_bgc ! MLD criterion
+      REAL(wp), INTENT(in   )                               :: k_maxchlinc    ! Max chl increment
+      LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n
+      LOGICAL,  INTENT(in   )                               :: ld_asmdin      ! Direct initialisation y/n
+      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)           :: chl_bkg        ! Surface chlorophyll
+      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) :: logchl_balinc  ! Balancing increments
+      !!
+      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
+      INTEGER                                               :: jkmax          ! Loop index
+      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
+      REAL(wp),                DIMENSION(jpi,jpj)           :: chl_inc        ! Chlorophyll increments
+      REAL(wp),                DIMENSION(jpi,jpj)           :: zmld           ! Mixed layer depth
+      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
+      !!---------------------------------------------------------------------------
+
+      ! 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 k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( chl_bkg(ji,jj) > 0.0 ) THEN
+               chl_inc(ji,jj) = 10**( LOG10( chl_bkg(ji,jj) ) + logchl_bkginc(ji,jj) ) - chl_bkg(ji,jj)
+               IF ( k_maxchlinc > 0.0 ) THEN
+                  chl_inc(ji,jj) = MAX( -1.0 * k_maxchlinc, MIN( chl_inc(ji,jj), k_maxchlinc ) )
+               ENDIF
+            ELSE
+               chl_inc(ji,jj) = 0.0
+            ENDIF
+         END DO
+      END DO
+      
+      ! Select mixed layer
+      IF ( ld_asmdin ) THEN
+         CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
+            &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
+         zmld(:,:) = mld_max_bkg(:,:)
+      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 LogChl assimilation,', &
+                  &           ' but ln_kara=.false.' )
+            ENDIF
+#else
+            CALL ctl_stop( ' Kara mixed layer requested for LogChl 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 LogChl 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
+      
+      IF ( ld_logchlbal ) 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(aincper), 1, kmt(:,:), tmask(:,:,:),                    &
+            &               zmld(:,:), mld_max_bkg(:,:), chl_inc(:,:), 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
+         logchl_balinc(:,:,:,jp_had_nut) = outincs(:,:,:,i_tracer(1))
+         logchl_balinc(:,:,:,jp_had_phy) = outincs(:,:,:,i_tracer(2))
+         logchl_balinc(:,:,:,jp_had_zoo) = outincs(:,:,:,i_tracer(3))
+         logchl_balinc(:,:,:,jp_had_pdn) = outincs(:,:,:,i_tracer(4))
+         logchl_balinc(:,:,:,jp_had_dic) = outincs(:,:,:,i_tracer(5))
+         logchl_balinc(:,:,:,jp_had_alk) = outincs(:,:,:,i_tracer(6))
+      
+      ELSE   ! No nitrogen balancing
+      
+         ! Initialise phytoplankton increment to zero
+         logchl_balinc(:,:,:,jp_had_phy) = 0.0
+         
+         ! Convert surface chlorophyll increment to phytoplankton nitrogen
+         logchl_balinc(:,:,1,jp_had_phy) = ( cchl_p_bkg(:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:)
+         
+         ! 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
+               !
+               DO jk = 2, jkmax
+                  logchl_balinc(ji,jj,jk,jp_had_phy) = logchl_balinc(ji,jj,1,jp_had_phy)
+               END DO
+               !
+            END DO
+         END DO
+
+         ! Set other balancing increments to zero
+         logchl_balinc(:,:,:,jp_had_nut) = 0.0
+         logchl_balinc(:,:,:,jp_had_zoo) = 0.0
+         logchl_balinc(:,:,:,jp_had_pdn) = 0.0
+         logchl_balinc(:,:,:,jp_had_dic) = 0.0
+         logchl_balinc(:,:,:,jp_had_alk) = 0.0
+      
+      ENDIF
+
+   END SUBROUTINE asm_logchl_bal_hadocc
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option : Empty routine
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, &
+      &                              k_maxchlinc, ld_logchlbal,              &
+      &                              pgrow_avg_bkg, ploss_avg_bkg,           &
+      &                              phyt_avg_bkg, mld_max_bkg,              &
+      &                              chl_bkg, cchl_p_bkg,                    &
+      &                              tracer_bkg, logchl_balinc )
+      REAL    :: logchl_bkginc(:,:)
+      REAL    :: aincper
+      INTEGER :: mld_choice_bgc
+      REAL    :: k_maxchlinc
+      LOGICAL :: ld_logchlbal
+      REAL    :: pgrow_avg_bkg(:,:)
+      REAL    :: ploss_avg_bkg(:,:)
+      REAL    :: phyt_avg_bkg(:,:)
+      REAL    :: mld_max_bkg(:,:)
+      REAL    :: chl_bkg(:,:)
+      REAL    :: cchl_p_bkg(:,:)
+      REAL    :: tracer_bkg(:,:,:,:)
+      REAL    :: logchl_balinc(:,:,:,:)
+      WRITE(*,*) 'asm_logchl_bal_hadocc: You should not have seen this print! error?'
+   END SUBROUTINE asm_logchl_bal_hadocc
+#endif
+
+   !!======================================================================
+END MODULE asmlogchlbal_hadocc
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90	(revision 9292)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90	(revision 9292)
@@ -0,0 +1,439 @@
+MODULE asmlogchlbal_medusa
+   !!======================================================================
+   !!                       ***  MODULE asmlogchlbal_medusa  ***
+   !! Calculate increments to MEDUSA based on surface logchl 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 asmlogchlbal_hadocc
+   !!----------------------------------------------------------------------
+#if defined key_asminc && defined key_medusa && defined key_foam_medusa
+   !!----------------------------------------------------------------------
+   !! 'key_asminc'          : assimilation increment interface
+   !! 'key_medusa'          : MEDUSA model
+   !! 'key_foam_medusa'     : MEDUSA extras for FOAM OBS and ASM
+   !!----------------------------------------------------------------------
+   !! asm_logchl_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 zdfmxl                              ! mixed layer depth
+   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_logchl_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_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, &
+      &                              k_maxchlinc, ld_logchlbal, ld_asmdin,   &
+      &                              pgrow_avg_bkg, ploss_avg_bkg,           &
+      &                              phyt_avg_bkg, mld_max_bkg,              &
+      &                              tracer_bkg, logchl_balinc )
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_logchl_bal_medusa  ***
+      !!
+      !! ** Purpose :   calculate increments to MEDUSA from logchl increments
+      !!
+      !! ** Method  :   convert logchl increments to chl increments
+      !!                average up MEDUSA to look like HadOCC
+      !!                call nitrogen balancing scheme
+      !!                separate back out to MEDUSA
+      !!
+      !! ** Action  :   populate logchl_balinc
+      !!
+      !! References :   Hemmings et al., 2008, J. Mar. Res.
+      !!                Ford et al., 2012, Ocean Sci.
+      !!---------------------------------------------------------------------------
+      !!
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: logchl_bkginc  ! logchl increments
+      REAL(wp), INTENT(in   )                               :: aincper        ! Assimilation period
+      INTEGER,  INTENT(in   )                               :: mld_choice_bgc ! MLD criterion
+      REAL(wp), INTENT(in   )                               :: k_maxchlinc    ! Max chl increment
+      LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n
+      LOGICAL,  INTENT(in   )                               :: ld_asmdin      ! Direct initialisation y/n
+      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) :: logchl_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_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_dtc_det   ! Ratio of jpdtc:jpdet
+      REAL(wp),                DIMENSION(jpi,jpj)           :: chl_inc        ! Chlorophyll increments
+      REAL(wp),                DIMENSION(jpi,jpj)           :: medusa_chl     ! MEDUSA total chlorophyll
+      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
+      REAL(wp),                DIMENSION(jpi,jpj)           :: zmld           ! Mixed layer depth
+      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
+      !!---------------------------------------------------------------------------
+
+      ! 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 k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      medusa_chl(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( medusa_chl(ji,jj) > 0.0 ) THEN
+               chl_inc(ji,jj) = 10**( LOG10( medusa_chl(ji,jj) ) + logchl_bkginc(ji,jj) ) - medusa_chl(ji,jj)
+               IF ( k_maxchlinc > 0.0 ) THEN
+                  chl_inc(ji,jj) = MAX( -1.0 * k_maxchlinc, MIN( chl_inc(ji,jj), k_maxchlinc ) )
+               ENDIF
+            ELSE
+               chl_inc(ji,jj) = 0.0
+            ENDIF
+         END DO
+      END DO
+      
+      ! Select mixed layer
+      IF ( ld_asmdin ) THEN
+         CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
+            &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
+         zmld(:,:) = mld_max_bkg(:,:)
+      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 LogChl assimilation,', &
+                  &           ' but ln_kara=.false.' )
+            ENDIF
+#else
+            CALL ctl_stop( ' Kara mixed layer requested for LogChl 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 LogChl 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
+      
+      IF ( ld_logchlbal ) 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(aincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
+            &               zmld(:,:), mld_max_bkg(:,:), chl_inc(:,:), 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
+         logchl_balinc(:,:,:,:) = 0.0
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) ) THEN
+                     ! Phytoplankton nitrogen and silicate 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 = 1.0 - zfrac_phn
+                     zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd)
+                     logchl_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn
+                     logchl_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd
+                     logchl_balinc(ji,jj,jk,jppds) = logchl_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
+
+                     ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen
+                     ! Not using chl_inc directly as it's only 2D
+                     ! This method should give same results at surface as splitting chl_inc 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)
+                     logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
+                     logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,jk,jpphd) * zrat_chd_phd
+                  ENDIF
+
+                  IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN
+                     ! Zooplankton nitrogen split up based on existing ratios
+                     zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
+                     zfrac_zme = 1.0 - zfrac_zmi
+                     logchl_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi
+                     logchl_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme
+                  ENDIF
+
+                  ! Nitrogen nutrient straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1))
+
+                  ! Nitrogen detritus straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4))
+
+                  ! DIC straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5))
+
+                  ! Alkalinity straight from balancing scheme
+                  logchl_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) - logchl_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
+                     logchl_balinc(ji,jj,jk,jpsil) = logchl_balinc(ji,jj,jk,jppds) * (-1.0)
+                  ENDIF
+
+                  IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN
+                     ! Carbon detritus based on existing ratios
+                     zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet)
+                     logchl_balinc(ji,jj,jk,jpdtc) = logchl_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
+                  ENDIF
+
+                  ! Do nothing with iron or oxygen for the time being
+                  logchl_balinc(ji,jj,jk,jpfer) = 0.0
+                  logchl_balinc(ji,jj,jk,jpoxy) = 0.0
+                  
+               END DO
+            END DO
+         END DO
+      
+      ELSE   ! No nitrogen balancing
+      
+         ! Initialise individual chlorophyll increments to zero
+         logchl_balinc(:,:,:,jpchn) = 0.0
+         logchl_balinc(:,:,:,jpchd) = 0.0
+         
+         ! Split up total surface chlorophyll increments
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( medusa_chl(ji,jj) > 0.0 ) THEN
+                  zfrac_chn = tracer_bkg(ji,jj,1,jpchn) / medusa_chl(ji,jj)
+                  zfrac_chd = 1.0 - zfrac_chn
+                  logchl_balinc(ji,jj,1,jpchn) = chl_inc(ji,jj) * zfrac_chn
+                  logchl_balinc(ji,jj,1,jpchd) = chl_inc(ji,jj) * zfrac_chd
+               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 ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                     jkmax = jk
+                  ENDIF
+               END DO
+               !
+               DO jk = 2, jkmax
+                  logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,1,jpchn)
+                  logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,1,jpchd)
+               END DO
+               !
+            END DO
+         END DO
+
+         ! Set other balancing increments to zero
+         logchl_balinc(:,:,:,jpphn) = 0.0
+         logchl_balinc(:,:,:,jpphd) = 0.0
+         logchl_balinc(:,:,:,jppds) = 0.0
+         logchl_balinc(:,:,:,jpzmi) = 0.0
+         logchl_balinc(:,:,:,jpzme) = 0.0
+         logchl_balinc(:,:,:,jpdin) = 0.0
+         logchl_balinc(:,:,:,jpsil) = 0.0
+         logchl_balinc(:,:,:,jpfer) = 0.0
+         logchl_balinc(:,:,:,jpdet) = 0.0
+         logchl_balinc(:,:,:,jpdtc) = 0.0
+         logchl_balinc(:,:,:,jpdic) = 0.0
+         logchl_balinc(:,:,:,jpalk) = 0.0
+         logchl_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
+               logchl_balinc(:,:,jk,jn) = logchl_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
+            END DO
+         END DO
+      ENDIF
+
+   END SUBROUTINE asm_logchl_bal_medusa
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option : Empty routine
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE asm_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, &
+      &                              k_maxchlinc, ld_logchlbal,              &
+      &                              pgrow_avg_bkg, ploss_avg_bkg,           &
+      &                              phyt_avg_bkg, mld_max_bkg,              &
+      &                              tracer_bkg, logchl_balinc )
+      REAL    :: logchl_bkginc(:,:)
+      REAL    :: aincper
+      INTEGER :: mld_choice_bgc
+      REAL    :: k_maxchlinc
+      LOGICAL :: ld_logchlbal
+      REAL    :: pgrow_avg_bkg(:,:)
+      REAL    :: ploss_avg_bkg(:,:)
+      REAL    :: phyt_avg_bkg(:,:)
+      REAL    :: mld_max_bkg(:,:)
+      REAL    :: tracer_bkg(:,:,:,:)
+      REAL    :: logchl_balinc(:,:,:,:)
+      WRITE(*,*) 'asm_logchl_bal_medusa: You should not have seen this print! error?'
+   END SUBROUTINE asm_logchl_bal_medusa
+#endif
+
+   !!======================================================================
+END MODULE asmlogchlbal_medusa
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90	(revision 9292)
@@ -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_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90	(revision 9292)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90	(revision 9292)
@@ -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_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 9292)
@@ -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_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 9292)
@@ -61,4 +61,5 @@
    USE asminc          ! assimilation increments     
    USE asmbkg          ! writing out state trajectory
+   USE asmbal          ! writing out biogeochemical assimilation balancing increments
    USE diaptr          ! poleward transports           (dia_ptr_init routine)
    USE diadct          ! sections transports           (dia_dct_init routine)
@@ -160,4 +161,6 @@
                 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics
                 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
+                IF( ln_logchlinc ) CALL logchl_asm_inc( nit000 - 1 )
+                IF( ln_fco2inc .OR. ln_pco2inc ) CALL pco2_asm_inc( nit000 - 1 )
              ENDIF
           ENDIF
@@ -191,4 +194,6 @@
 
       IF( lk_diaobs   )   CALL dia_obs_wri
+      !
+      IF( ( lk_asminc ).AND.( ln_balwri ) ) CALL asm_bal_wri( nitend )  ! Output balancing increments
       !
       IF( ln_icebergs )   CALL icb_end( nitend )
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 9292)
@@ -253,4 +253,8 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( lk_asminc .AND. ln_asmiau .AND. ln_logchlinc ) &
+         &               CALL logchl_asm_inc( kstp )  ! logchl assimilation
+      IF( lk_asminc .AND. ln_asmiau .AND. ( ln_fco2inc .OR. ln_pco2inc ) ) &
+         &               CALL pco2_asm_inc( kstp )    ! pCO2/fCO2 assimilation
                          CALL trc_stp( kstp )         ! time-stepping
 #endif
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/air_sea.F90	(revision 9292)
@@ -77,4 +77,7 @@
                                    zn_dms_chd, zn_dms_chn, zn_dms_din,    &
                                    zn_dms_mld, zn_dms_qsr,                &
+# if defined key_foam_medusa
+                                   f2_pco2w, f2_fco2w,                    &
+# endif
                                    xnln, xnld 
       USE trc,               ONLY: med_diag
@@ -87,4 +90,7 @@
 #  else
       USE trcco2_medusa,     ONLY: trc_co2_medusa
+#   if defined key_foam_medusa
+      USE mocsy_mainmod,     ONLY: p2fCO2
+#   endif
 #  endif
       USE trcdms_medusa,     ONLY: trc_dms_medusa
@@ -330,4 +336,8 @@
                      iters, ' AT (', ji, ', ', jj, ', 1) AT ', kt
                endif
+#    if defined key_foam_medusa
+               !! DAF (Aug 2017): calculate fCO2 for observation operator
+               CALL p2fCO2( f_pco2w, ztmp, f_pp0, 0.0, 1, f_fco2w )
+#    endif
             ENDIF
          ENDDO
@@ -509,4 +519,9 @@
                               CO2flux_conv
                !! ENDIF
+#   if defined key_foam_medusa
+               !! 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_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 9292)
@@ -48,4 +48,8 @@
                                    f3_co3, f3_h2co3, f3_hco3,           &
                                    f3_omarg, f3_omcal, f3_pH,           &
+# if defined key_foam_medusa
+                                   mld_max, pgrow_avg,                  &
+                                   ploss_avg, phyt_avg,                 &
+# endif
                                    za_sed_c, za_sed_ca, za_sed_fe,      &
                                    za_sed_n, za_sed_si,                 &
@@ -57,4 +61,7 @@
       USE trc,               ONLY: med_diag, nittrc000, trn 
       USE trcnam_trp,        ONLY: ln_trcadv_cen2, ln_trcadv_tvd
+# if defined key_foam_medusa
+      USE zdfmxl,            ONLY: hmld
+# endif
  
       !! time (integer timestep)
@@ -115,4 +122,25 @@
 # endif      
 
+# if defined key_foam_medusa
+!!----------------------------------------------------------------------
+!! 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_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 9292)
@@ -35,5 +35,9 @@
       USE bio_medusa_mod
       USE par_oce,           ONLY: jpi, jpj, jpk
+# if defined key_foam_medusa
+      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max
+# else
       USE sms_medusa,        ONLY: jdms
+# endif
       USE trc,               ONLY: ln_diatrc, med_diag, nittrc000 
       USE in_out_manager,    ONLY: lwp
@@ -192,4 +196,13 @@
       fslowsinkc(:,:) = 0.0
 # endif      
+      !!
+# if defined key_foam_medusa
+      pgrow_avg(:,:) = 0.0
+      ploss_avg(:,:) = 0.0
+      phyt_avg(:,:)  = 0.0
+      IF( kt == nittrc000 ) THEN
+         mld_max(:,:) = 0.0
+      ENDIF
+# endif
       !!
       !! allocate and initiate 2D diag
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 9292)
@@ -36,13 +36,33 @@
                                    fdpn, fdpn2, fdzme, fdzme2,             &
                                    fdzmi, fdzmi2, fsdiss, fsin,            &
+# if defined key_foam_medusa
+                                   fdep1, fprn, fprd,                      &
+                                   fgmepd, fgmepn, fgmipn,                 &
+# endif
                                    zphd, zphn, zpds, zzme, zzmi
+# if defined key_foam_medusa
+      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask
+      USE par_kind,          ONLY: wp
+# else
       USE dom_oce,           ONLY: tmask
+# endif
       USE par_oce,           ONLY: jpim1, jpjm1
       USE phytoplankton_mod, ONLY: phytoplankton
       USE sms_medusa,        ONLY: jmpd, jmpn, jmzme, jmzmi,               &
+# if defined key_foam_medusa
+                                   pgrow_avg, ploss_avg, phyt_avg,         &
+# endif
                                    xkphd, xkphn, xkzme, xkzmi,             &
                                    xmetapd, xmetapn, xmetazme, xmetazmi,   &
                                    xmpd, xmpn, xmzme, xmzmi, xsdiss
+# if defined key_foam_medusa
+      USE zdfmxl,            ONLY: hmld
+# endif
       USE zooplankton_mod,   ONLY: zooplankton
+
+# if defined key_foam_medusa
+   !!* Substitution
+#  include "domzgr_substitute.h90"
+# endif
 
       !! Level
@@ -50,4 +70,8 @@
 
       INTEGER :: ji, jj
+
+# if defined key_foam_medusa
+      REAL(wp) :: fq0
+# endif
 
       !!-------------------------------------------------------------------
@@ -188,4 +212,39 @@
       ENDDO
 
+# if defined key_foam_medusa
+      !! 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_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 9292)
@@ -212,4 +212,10 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: f2_ccd_arg  !: 2D aragonite CCD depth
 !!
+#if defined key_foam_medusa
+!! 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
+!!
+#endif
 !! 2D fields of organic and inorganic material sedimented on the seafloor
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: zb_sed_n    !: 2D organic nitrogen   (before)
@@ -391,4 +397,7 @@
       & 926.4960, 935.7040 /
 # endif
+# if defined key_foam_medusa
+   REAL(wp) ::   xobs_xco2a   !: Observed atmospheric xCO2, read in
+# endif
 #endif
 
@@ -441,4 +450,16 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: cmask       !: ???
 
+#if defined key_foam_medusa
+!!----------------------------------------------------------------------
+!! Parameters required for ocean colour assimilation
+!!----------------------------------------------------------------------
+!!
+   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
+!!
+#endif
+
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -453,5 +474,5 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_warn
-      INTEGER ::   ierr(8)        ! Local variables
+      INTEGER ::   ierr(9)        ! Local variables
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -463,4 +484,7 @@
       !* 2D and 3D fields of carbonate system parameters
       ALLOCATE( f2_ccd_cal(jpi,jpj)  , f2_ccd_arg(jpi,jpj)  ,       &
+#  if defined key_foam_medusa
+                f2_pco2w(jpi,jpj)    , f2_fco2w(jpi,jpj)    ,       &
+#  endif
          &      f3_pH(jpi,jpj,jpk)   , f3_h2co3(jpi,jpj,jpk),       &
          &      f3_hco3(jpi,jpj,jpk) , f3_co3(jpi,jpj,jpk)  ,       &
@@ -511,4 +535,8 @@
          &      ffln(jpi,jpj,jpk)    , fflf(jpi,jpj,jpk)    ,       &
          &      ffls(jpi,jpj,jpk)    , cmask(jpi,jpj)       ,    STAT=ierr(8) ) 
+# if defined key_foam_medusa
+      ALLOCATE( pgrow_avg(jpi,jpj)   , ploss_avg(jpi,jpj)   ,       &
+         &      phyt_avg(jpi,jpj)    , mld_max(jpi,jpj)     ,    STAT=ierr(9) )
+# endif
 #endif
       !
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 9292)
@@ -102,5 +102,11 @@
       !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
       USE sbc_oce,                    ONLY: lk_oasis
+# if defined key_foam_medusa
+      USE sms_medusa,                 ONLY: hist_pco2, xobs_xco2a,          &
+                                            pgrow_avg, ploss_avg,           &
+                                            phyt_avg, mld_max
+# else
       USE sms_medusa,                 ONLY: hist_pco2
+# endif
       USE trc,                        ONLY: ln_rsttr, nittrc000, trn
       USE bio_medusa_init_mod,        ONLY: bio_medusa_init
@@ -319,4 +325,12 @@
          f_xco2a(:,:) = fq4
       endif
+#  if defined key_foam_medusa
+      IF ( xobs_xco2a > 0.0 ) THEN
+         IF(lwp) WRITE(numout,*) ' using observed atm pCO2 = ', xobs_xco2a
+         f_xco2a(:,:) = xobs_xco2a
+      ELSE
+         IF(lwp) WRITE(numout,*) ' xobs_xco2a <= 0 so using default atm pCO2'
+      ENDIF
+#  endif
 #  if defined key_axy_pi_co2
       !! OCMIP pre-industrial pCO2
@@ -358,7 +372,13 @@
       !!          x * 30d + 1*rdt(i.e: mod = rdt)   
       !!          ++ need to pass carb-chem output var through restarts
+#if defined key_foam_medusa
+      !! DAF (Aug 2017): For FOAM we want to run daily
+      If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
+           (mod(kt*rdt,86400.) == rdt) ) THEN
+#else
       If ( (kt == nittrc000 .AND. .NOT.ln_rsttr) .OR.                        &
            ((86400*mod(nn_date0,100) + mod(kt*rdt,2592000.)) == rdt) ) THEN
-         !!---------------------------------------------------------------
+#endif
+         !!----------------------------------------------------------------------
          !! Calculate the carbonate chemistry for the whole ocean on the first
          !! simulation timestep and every month subsequently; the resulting 3D
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcnam_medusa.F90	(revision 9292)
@@ -90,4 +90,7 @@
 #if defined key_roam
       NAMELIST/natroam/ xthetaphy,xthetazoo,xthetanit,        &
+# if defined key_foam_medusa
+      &    xobs_xco2a,                                        &
+# endif
       &    xthetarem,xo2min 
 #endif
@@ -490,4 +493,11 @@
          &   ' key_debug_medusa                                                       = INACTIVE'
 # endif
+#if defined key_foam_medusa
+         WRITE(numout,*)     &
+         &   ' key_foam_medusa                                                        = ACTIVE'
+#else
+         WRITE(numout,*)     &
+         &   ' key_foam_medusa                                                        = INACTIVE'
+#endif
          WRITE(numout,*) ' '
 
@@ -1030,4 +1040,7 @@
       xthetarem = 0.
       xo2min    = 0.
+# if defined key_foam_medusa
+      xobs_xco2a = 0.
+# endif
 
       !READ(numnatm,natroam)
@@ -1049,4 +1062,7 @@
 !!       xthetarem :  oxygen consumption by carbon remineralisation
 !!       xo2min    :  oxygen minimum concentration
+# if defined key_foam_medusa
+!!       xobs_xco2a : observed atmospheric xCO2
+# endif
 
       IF(lwp) THEN
@@ -1066,4 +1082,8 @@
           WRITE(numout,*)     &
           &   ' oxygen minimum concentration                               xo2min      = ', xo2min
+# if defined key_foam_medusa
+          WRITE(numout,*)     &
+          &   ' observed atmospheric xCO2                                  xobs_xco2a  = ', xobs_xco2a
+# endif
        ENDIF
 
Index: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 9291)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 9292)
@@ -45,4 +45,7 @@
    USE oce,     ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl  !! Coupling variable
    USE trcstat
+#if defined key_foam_medusa
+   USE obs_const, ONLY: obfillflt  ! Observation operator fill value
+#endif
 
    IMPLICIT NONE
@@ -338,6 +341,34 @@
          IF(lwp) WRITE(numout,*) 'Or don t start from uncomplete restart...' 
       ENDIF
+      !
+#  if defined key_foam_medusa
+      !! 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 defined key_foam_medusa
+      !! 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,4 +541,18 @@
       call trc_rst_dia_stat( f2_ccd_arg(:,:),'CCD_ARG')
       !!
+#  if defined key_foam_medusa
+      !! 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 defined key_foam_medusa
+      !! 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
 !!
