Index: /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8427)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8428)
@@ -1293,4 +1293,6 @@
     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_temnofreeze = .false. !  Logical to not add increments if temperature would fall below freezing
     nitbkg    = 0          !  Timestep of background in [0,nitend-nit000-1]
@@ -1302,4 +1304,8 @@
     salfixmin = -9999      !  Minimum salinity after applying the increments
     nn_divdmp = 0          !  Number of iterations of divergence damping operator
+    mld_choice_bgc = 5     !  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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 8428)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 8428)
@@ -0,0 +1,116 @@
+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 && defined key_foam_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
+
+         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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8427)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8428)
@@ -22,4 +22,5 @@
    !!   ssh_asm_inc    : Apply the SSH increment
    !!   seaice_asm_inc : Apply the seaice increment
+   !!   logchl_asm_inc : Apply the logchl increment
    !!----------------------------------------------------------------------
    USE wrk_nemo         ! Memory Allocation
@@ -44,4 +45,29 @@
 #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 par_medusa
+#elif defined key_hadocc
+   USE asmlogchlbal_hadocc, ONLY: &
+      & asm_logchl_bal_hadocc
+   USE par_hadocc
+#endif
 
    IMPLICIT NONE
@@ -54,4 +80,5 @@
    PUBLIC   ssh_asm_inc    !: Apply the SSH increment
    PUBLIC   seaice_asm_inc !: Apply the seaice increment
+   PUBLIC   logchl_asm_inc !: Apply the logchl increment
 
 #if defined key_asminc
@@ -61,4 +88,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 +95,6 @@
    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_salfix = .FALSE.      !: Apply minimum salinity check
    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
@@ -91,4 +121,19 @@
    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
+#if defined key_top
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim
+#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    = 5   !: 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 +178,12 @@
       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_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 +210,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 +216,6 @@
          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 Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
          WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
@@ -176,4 +226,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 +275,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 ) )) &
+         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' and ln_logchlinc is set to .true.', &
          &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
          &                ' Inconsistent options')
@@ -233,6 +287,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 ) )  &
+         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' and ln_logchlinc are set to .false. :', &
          &                ' The assimilation increments are not applied')
 
@@ -254,4 +309,16 @@
          &                ' Background time step for Direct Initialization is outside', &
          &                ' the cycle interval')
+
+      IF ( ( ln_balwri ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN
+         CALL ctl_warn( ' Balancing increments are only calculated for logchl', &
+            &           ' Not assimilating logchl, 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 ( nstop > 0 ) RETURN       ! if there are any errors then go no further
@@ -350,5 +417,14 @@
       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_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
+         &  .OR.( ln_logchlinc ) ) THEN
 
          !--------------------------------------------------------------------
@@ -422,4 +498,13 @@
             ! to allow for differences in masks
             WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
+         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
 
@@ -1155,4 +1240,119 @@
 
    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()
+         CALL ctl_stop( 'Attempting to assimilate logchl into MEDUSA, ', &
+            &           'but not fully implemented yet' )
+#elif defined key_hadocc
+         CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, &
+            &                        rn_maxchlinc, ln_logchlbal, 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, as not a tendency
+#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
    
    !!======================================================================
Index: /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_hadocc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_hadocc.F90	(revision 8428)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_hadocc.F90	(revision 8428)
@@ -0,0 +1,293 @@
+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 trc,           ONLY: trn, trb,    & ! HadOCC variables
+      &                     HADOCC_CHL,  &
+      &                     pgrow_avg,   &
+      &                     ploss_avg,   &
+      &                     phyt_avg,    &
+      &                     mld_max
+   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, 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
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc  ! Balancing increments
+      !!
+      INTEGER                                               :: ji, jj, jk     ! 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 ( HADOCC_CHL(ji,jj,1) > 0.0 ) THEN
+               chl_inc(ji,jj) = 10**( LOG10( HADOCC_CHL(ji,jj,1) ) + logchl_bkginc(ji,jj) ) - HADOCC_CHL(ji,jj,1)
+               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
+      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
+      
+      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)) = trb(:,:,:,jp_had_nut)
+         bstate(:,:,:,i_tracer(2)) = trb(:,:,:,jp_had_phy)
+         bstate(:,:,:,i_tracer(3)) = trb(:,:,:,jp_had_zoo)
+         bstate(:,:,:,i_tracer(4)) = trb(:,:,:,jp_had_pdn)
+         bstate(:,:,:,i_tracer(5)) = trb(:,:,:,jp_had_dic)
+         bstate(:,:,:,i_tracer(6)) = trb(:,:,:,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(:,:), chl_inc(:,:), cchl_p(:,:,1),  &
+            &               nbal_active, phyt_avg(:,:),                            &
+            &               gl_active, pgrow_avg(:,:), ploss_avg(:,:),             &
+            &               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(:,:,1) / (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, logchl_balinc )
+      REAL    :: logchl_bkginc(:,:)
+      REAL    :: aincper
+      INTEGER :: mld_choice_bgc
+      REAL    :: k_maxchlinc
+      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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90	(revision 8427)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90	(revision 8428)
@@ -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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8427)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8428)
@@ -61,4 +61,5 @@
    USE asminc          ! assimilation increments     
    USE asmbkg          ! writing out state trajectory
+   USE asmbal          ! writing out assimilation balancing increments
    USE diaptr          ! poleward transports           (dia_ptr_init routine)
    USE diadct          ! sections transports           (dia_dct_init routine)
@@ -160,4 +161,5 @@
                 IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics
                 IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
+                IF( ln_logchlinc ) CALL logchl_asm_inc( nit000 - 1 )
              ENDIF
           ENDIF
@@ -191,4 +193,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_surf_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8427)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8428)
@@ -253,4 +253,6 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( lk_asminc .AND. ln_asmiau .AND. ln_logchlinc ) &
+         &               CALL logchl_asm_inc( kstp )  ! logchl assimilation
                          CALL trc_stp( kstp )         ! time-stepping
 #endif
