Index: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90
===================================================================
--- branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 9318)
+++ branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmbal.F90	(revision 9331)
@@ -77,8 +77,37 @@
          IF ( ln_logchltotinc .OR. ln_logchlpftinc ) THEN
 #if defined key_fabm
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chl1', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl1) )
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chl2', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl2) )
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chl3', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl3) )
-            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chl4', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl4) )
+            CALL iom_rstput( kt, kt, inum, 'bal_chl1', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl1) )
+            CALL iom_rstput( kt, kt, inum, 'bal_chl2', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl2) )
+            CALL iom_rstput( kt, kt, inum, 'bal_chl3', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl3) )
+            CALL iom_rstput( kt, kt, inum, 'bal_chl4', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl4) )
+            IF ( ln_logchlbal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'bal_p1n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p1c', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1c) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p1p', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1p) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p1s', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1s) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p2n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p2c', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2c) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p2p', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2p) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p3n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p3c', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3c) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p3p', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3p) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p4n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p4c', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4c) )
+               CALL iom_rstput( kt, kt, inum, 'bal_p4p', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4p) )
+               CALL iom_rstput( kt, kt, inum, 'bal_n3n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_n3n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_n4n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_n4n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_z4c', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_z4c) )
+               CALL iom_rstput( kt, kt, inum, 'bal_z5n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_z5n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_z5c', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_z5c) )
+               CALL iom_rstput( kt, kt, inum, 'bal_z5p', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_z5p) )
+               CALL iom_rstput( kt, kt, inum, 'bal_z6n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_z6n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_z6c', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_z6c) )
+               CALL iom_rstput( kt, kt, inum, 'bal_z6p', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_z6p) )
+               CALL iom_rstput( kt, kt, inum, 'bal_r4n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_r4n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_r6n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_r6n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_r8n', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_r8n) )
+               CALL iom_rstput( kt, kt, inum, 'bal_o3c', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3c) )
+               CALL iom_rstput( kt, kt, inum, 'bal_o3a', logchl_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3a) )
+            ENDIF
 #elif defined key_medusa && defined key_foam_medusa
             CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', logchl_balinc(:,:,:,jpchn) )
Index: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 9318)
+++ branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 9331)
@@ -51,4 +51,11 @@
 #endif
    USE asminc, ONLY: ln_avgbkg
+#if defined key_fabm
+   USE trc, ONLY:           & ! HadOCC diagnostic variables
+      & pgrow_avg,          &
+      & ploss_avg,          &
+      & phyt_avg,           &
+      & mld_max
+#endif
    IMPLICIT NONE
    PRIVATE
@@ -196,4 +203,10 @@
 #endif
                CALL iom_rstput( kt, nitavgbkg_r, inum, 'avt'    , avt_tavg)
+#if defined key_fabm
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'pgrow_avg'   , pgrow_avg             )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'ploss_avg'   , ploss_avg             )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'phyt_avg'    , phyt_avg              )
+               CALL iom_rstput( kt, nitavgbkg_r, inum, 'mld_max'     , mld_max               )
+#endif
                !
             ELSE
@@ -219,4 +232,10 @@
 #endif
                CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )
+#if defined key_fabm
+               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               )
+#endif
                !
             ENDIF
Index: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 9318)
+++ branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 9331)
@@ -101,4 +101,5 @@
    LOGICAL, PUBLIC :: ln_logchltotinc = .FALSE. !: No total log10(chlorophyll) increment
    LOGICAL, PUBLIC :: ln_logchlpftinc = .FALSE. !: No PFT   log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_logchlbal = .FALSE.
    LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
@@ -131,4 +132,8 @@
    REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim
 #endif
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg   ! Background phyto growth
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg   ! Background phyto loss
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg    ! Background phyto conc
+   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg     ! Background max MLD
    REAL(wp) :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation
                                       !: <= 0 implies no maximum applied (switch turned off)
@@ -206,5 +211,5 @@
       NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri, ln_avgbkg,                &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
-         &                 ln_logchltotinc, ln_logchlpftinc,               &
+         &                 ln_logchltotinc, ln_logchlpftinc, ln_logchlbal, &
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
@@ -279,4 +284,5 @@
          WRITE(numout,*) '      Logical switch for applying total logchl incs      ln_logchltotinc = ', ln_logchltotinc
          WRITE(numout,*) '      Logical switch for applying PFT   logchl incs      ln_logchlpftinc = ', ln_logchlpftinc
+         WRITE(numout,*) '      Logical switch for applying balancing logchl incs     ln_logchlbal = ', ln_logchlbal
          WRITE(numout,*) '      Number of logchl PFTs assimilated                       nn_asmpfts = ', nn_asmpfts
          WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
@@ -769,4 +775,30 @@
 
       ENDIF
+
+      IF ( ln_logchltotinc .OR. ln_logchlpftinc ) THEN
+
+         CALL iom_open( c_asmbkg, inum )
+
+         ALLOCATE( pgrow_avg_bkg(jpi,jpj)        )
+         ALLOCATE( ploss_avg_bkg(jpi,jpj)        )
+         ALLOCATE( phyt_avg_bkg(jpi,jpj)         )
+         ALLOCATE( mld_max_bkg(jpi,jpj)          )
+         pgrow_avg_bkg(:,:)  = 0.0
+         ploss_avg_bkg(:,:)  = 0.0
+         phyt_avg_bkg(:,:)   = 0.0
+         mld_max_bkg(:,:)    = 0.0
+
+         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)
+
+         CALL iom_close( inum )
+
+      ENDIF
       !
    END SUBROUTINE asm_inc_init
@@ -1412,5 +1444,8 @@
 #if defined key_fabm
          CALL asm_logchl_bal_ersem( ln_logchlpftinc, nn_asmpfts, mld_choice_bgc, &
-            &                       rn_maxchlinc, logchl_bkginc, logchl_balinc )
+            &                       rn_maxchlinc, logchl_bkginc, logchl_balinc,   &
+            &                       pgrow_avg_bkg, ploss_avg_bkg,           &
+            &                       phyt_avg_bkg, mld_max_bkg,              &
+            &                       ln_logchlbal, zincper )
 #elif defined key_medusa && defined key_foam_medusa
          !CALL asm_logchl_bal_medusa()
Index: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90
===================================================================
--- branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90	(revision 9318)
+++ branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90	(revision 9331)
@@ -21,4 +21,5 @@
    USE par_fabm                          ! FABM parameters
    USE par_trc,     ONLY: jptra          ! Tracer parameters
+   USE bioanalysis
 
    IMPLICIT NONE
@@ -27,8 +28,42 @@
    PUBLIC asm_logchl_bal_ersem
 
+   ! Default values for biological assimilation parameters
+   ! Should match Hemmings et al. (2008)
+   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
+   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
+   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
+   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
+   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
+   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
+   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
+   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
+   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
+   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
+   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
+   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
+   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
+   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
+   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
+   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
+   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
+   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
+   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
+   !
+   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
+   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
+   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
+   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
+   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
+   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
+   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
+   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
+
 CONTAINS
 
    SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
-      &                             k_maxchlinc, logchl_bkginc, logchl_balinc )
+      &                             k_maxchlinc, logchl_bkginc, logchl_balinc,   &
+      &                             pgrow_avg_bkg, ploss_avg_bkg,           &
+      &                             phyt_avg_bkg, mld_max_bkg,              &
+      &                             ld_logchlbal, aincper )
       !!---------------------------------------------------------------------------
       !!                    ***  ROUTINE asm_logchl_bal_ersem  ***
@@ -52,4 +87,10 @@
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts)     :: logchl_bkginc
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc
+      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
+      LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n
+      REAL(wp), INTENT(in   )                               :: aincper
       !!
       INTEGER                      :: ji, jj, jk
@@ -57,4 +98,42 @@
       REAL(wp)                     :: chl_tot, chl_inc
       REAL(wp), DIMENSION(jpi,jpj) :: zmld
+      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
+      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
+      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
+      REAL(wp),                DIMENSION(jpi,jpj)           :: chlinctot
+      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_n3n      ! Fraction of jp_fabm_n3n
+      REAL(wp)                                              :: zfrac_n4n      ! Fraction of jp_fabm_n4n
+      REAL(wp)                                              :: zfrac_r4n      ! Fraction of jp_fabm_r4n
+      REAL(wp)                                              :: zfrac_r6n      ! Fraction of jp_fabm_r6n
+      REAL(wp)                                              :: zfrac_r8n      ! Fraction of jp_fabm_r8n
+      REAL(wp)                                              :: zfrac_z4n      ! Fraction of jp_fabm_z4n
+      REAL(wp)                                              :: zfrac_z5n      ! Fraction of jp_fabm_z5n
+      REAL(wp)                                              :: zfrac_z6n      ! Fraction of jp_fabm_z6n
+      REAL(wp)                                              :: zfrac_p1n      ! Fraction of jp_fabm_p1n
+      REAL(wp)                                              :: zfrac_p2n      ! Fraction of jp_fabm_p2n
+      REAL(wp)                                              :: zfrac_p3n      ! Fraction of jp_fabm_p3n
+      REAL(wp)                                              :: zfrac_p4n      ! Fraction of jp_fabm_p4n
+      REAL(wp)                                              :: zrat_z4c_z4n
+      REAL(wp)                                              :: zrat_z5c_z5n
+      REAL(wp)                                              :: zrat_z5p_z5n
+      REAL(wp)                                              :: zrat_z6c_z6n
+      REAL(wp)                                              :: zrat_z6p_z6n
+      REAL(wp)                                              :: zrat_p1c_p1n
+      REAL(wp)                                              :: zrat_p1p_p1n
+      REAL(wp)                                              :: zrat_p1s_p1n
+      REAL(wp)                                              :: zrat_p2c_p2n
+      REAL(wp)                                              :: zrat_p2p_p2n
+      REAL(wp)                                              :: zrat_p3c_p3n
+      REAL(wp)                                              :: zrat_p3p_p3n
+      REAL(wp)                                              :: zrat_p4c_p4n
+      REAL(wp)                                              :: zrat_p4p_p4n
       !!---------------------------------------------------------------------------
 
@@ -176,6 +255,275 @@
          zmld(:,:) = hmlpt(:,:)
       END SELECT
+
+      IF ( ld_logchlbal ) THEN
+      
+         ! Set up model parameters to be passed into Hemmings balancing routine.
+         ! For now these are hardwired to the standard HadOCC parameter values
+         ! as this is what the scheme was developed for.
+         ! Obviously, HadOCC and ERSEM are rather different models, so this
+         ! isn't ideal, but there's not always direct analogues between the two
+         ! parameter sets, so it's the easiest way to get something running.
+         ! In the longer term, some serious MarMOT-based development is required.
+         modparm(1)  = 0.1                               ! grow_sat
+         modparm(2)  = 2.0                               ! psmax
+         modparm(3)  = 0.845                             ! par
+         modparm(4)  = 0.02                              ! alpha
+         modparm(5)  = 0.05                              ! resp_rate
+         modparm(6)  = 0.05                              ! pmort_rate
+         modparm(7)  = 0.01                              ! phyto_min
+         modparm(8)  = 0.05                              ! z_mort_1
+         modparm(9)  = 1.0                               ! z_mort_2
+         modparm(10) = 6.625                             ! c2n_p
+         modparm(11) = 5.625                             ! c2n_z
+         modparm(12) = 7.5                               ! c2n_d
+         modparm(13) = 0.01                              ! graze_threshold
+         modparm(14) = 2.0                               ! holling_coef
+         modparm(15) = 0.5                               ! graze_sat
+         modparm(16) = 2.0                               ! graze_max
+
+         ! Set up assimilation parameters to be passed into balancing routine
+         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
+         assimparm(2)  = balnutext
+         assimparm(3)  = balnutmin
+         assimparm(4)  = r
+         assimparm(5)  = beta_g
+         assimparm(6)  = beta_l
+         assimparm(7)  = beta_m
+         assimparm(8)  = a_g
+         assimparm(9)  = a_l
+         assimparm(10) = a_m
+         assimparm(11) = zfracb0
+         assimparm(12) = zfracb1
+         assimparm(13) = qrfmax
+         assimparm(14) = qafmax
+         assimparm(15) = zrfmax
+         assimparm(16) = zafmax
+         assimparm(17) = prfmax
+         assimparm(18) = incphymin
+         assimparm(19) = integnstep
+         assimparm(20) = pthreshold
+
+         ! Set up external tracer indices array bstate
+         i_tracer(1) = 1   ! nutrient
+         i_tracer(2) = 2   ! phytoplankton
+         i_tracer(3) = 3   ! zooplankton
+         i_tracer(4) = 4   ! detritus
+         i_tracer(5) = 5   ! DIC
+         i_tracer(6) = 6   ! Alkalinity
+
+         ! Set background state
+         bstate(:,:,:,i_tracer(1)) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + trn(:,:,:,jp_fabm_m1+jp_fabm_n4n)
+         bstate(:,:,:,i_tracer(2)) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1n) + trn(:,:,:,jp_fabm_m1+jp_fabm_p2n) + &
+            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_p3n) + trn(:,:,:,jp_fabm_m1+jp_fabm_p4n)
+         ! Z4c needs converting by qnc, hardwire for now
+         bstate(:,:,:,i_tracer(3)) = (trn(:,:,:,jp_fabm_m1+jp_fabm_z4c) * 0.0126 ) + &
+            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_z5n) + trn(:,:,:,jp_fabm_m1+jp_fabm_z6n)
+         bstate(:,:,:,i_tracer(4)) = trn(:,:,:,jp_fabm_m1+jp_fabm_r4n) + trn(:,:,:,jp_fabm_m1+jp_fabm_r6n) + &
+            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_r8n)
+         bstate(:,:,:,i_tracer(5)) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)
+         bstate(:,:,:,i_tracer(6)) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3a)
+
+         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
+         ! and nitrogen to biomass equivalent for PZD
+         ! Need a single number, so base on HadOCC
+         cchl_p(:,:) = 0.0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
+                  &   trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) .GT. 0.0 ) THEN
+                  cchl_p(ji,jj) = ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_p1c)  + trn(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) + &
+                     &              trn(ji,jj,1,jp_fabm_m1+jp_fabm_p3c)  + trn(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) ) /  &
+                     &            ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
+                     &              trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
+               ENDIF
+            END DO
+         END DO
+         n2be_p = ( 14.01 + ( 12.01 * 6.625 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
+         n2be_z = ( 14.01 + ( 12.01 * 5.625 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
+         n2be_d = ( 14.01 + ( 12.01 * 7.5 )   ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
+
+         chlinctot(:,:) = logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl1) + logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
+            &             logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl3) + logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl4)
+
+         ! 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(:,:), chlinctot(:,:), cchl_p(:,:), &
+            &               nbal_active, phyt_avg_bkg(:,:),                         &
+            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      &
+            &               subsurf_active, deepneg_active,                         &
+            &               deeppos_active, nutprof_active,                         &
+            &               bstate, outincs,                                        &
+            &               diag_active, diag,                                      &
+            &               diag_fulldepth_active, diag_fulldepth )
+         
+         ! Loop over each grid point partioning the increments
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  ! Nitrogen phytoplankton from balancing scheme
+                  ! Split according to current ratios [ChlTot] or assimilation [PFTs]
+                  ! Update carbon, phosphorus and silicon according to current ratios
+                  ! Already have chlorophyll
+                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) > 0.0 ) .AND. &
+                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) > 0.0 ) .AND. &
+                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) > 0.0 ) .AND. &
+                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) > 0.0 ) ) THEN
+
+                     IF ( ld_logchlpftinc ) THEN
+                        IF ( (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
+                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
+                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
+                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)) > 0.0 ) THEN
+                           zfrac_p1n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / &
+                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
+                           zfrac_p2n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / &
+                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
+                           zfrac_p3n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / &
+                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
+                           zfrac_p4n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / &
+                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
+                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
+                        ELSE
+                           zfrac_p1n = 0.25
+                           zfrac_p2n = 0.25
+                           zfrac_p3n = 0.25
+                           zfrac_p4n = 0.25
+                        ENDIF
+                     ELSE
+                        zfrac_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) / &
+                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
+                        zfrac_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) / &
+                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
+                        zfrac_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) / &
+                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
+                        zfrac_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) / &
+                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
+                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
+                     ENDIF
+                     
+                     zrat_p1c_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
+                     zrat_p1p_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
+                     zrat_p1s_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
+                     zrat_p2c_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
+                     zrat_p2p_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
+                     zrat_p3c_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
+                     zrat_p3p_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
+                     zrat_p4c_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
+                     zrat_p4p_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
+                     
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1c_p1n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1p_p1n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1s_p1n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2c_p2n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2p_p2n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3c_p3n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3p_p3n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4c_p4n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4p_p4n
+
+                  ENDIF
+
+                  ! Nitrogen nutrient from balancing scheme
+                  ! Split between nitrate and ammonium according to current ratios
+                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) > 0.0 ) .AND. &
+                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) > 0.0 ) ) THEN
+                     zfrac_n3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) / &
+                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) + trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n))
+                     zfrac_n4n = 1.0 - zfrac_n3n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n
+                  ENDIF
+
+                  ! Nitrogen zooplankton from balancing scheme
+                  ! Split according to current ratios
+                  ! Update carbon and phosphorus according to current ratios
+                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) > 0.0 ) .AND. &
+                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) > 0.0 ) .AND. &
+                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) > 0.0 ) ) THEN
+                     zfrac_z4n = (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) / &
+                        &        ((trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) + &
+                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
+                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n))
+                     zfrac_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) / &
+                        &        ((trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) + &
+                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
+                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n))
+                     zfrac_z6n = 1.0 - zfrac_z4n - zfrac_z5n
+                     zrat_z4c_z4n = 1.0 / 0.0126
+                     zrat_z5c_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n)
+                     zrat_z5p_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n)
+                     zrat_z6c_z6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)
+                     zrat_z6p_z6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5c_z5n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6c_z6n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5p_z5n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6p_z6n
+                  ENDIF
+
+                  ! Nitrogen detritus from balancing scheme
+                  ! Split according to current ratios
+                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) > 0.0 ) .AND. &
+                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) > 0.0 ) .AND. &
+                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) > 0.0 ) ) THEN
+                     zfrac_r4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) / &
+                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
+                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
+                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
+                     zfrac_r6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) / &
+                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
+                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
+                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
+                     zfrac_r8n = 1.0 - zfrac_r4n - zfrac_r6n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n
+                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n
+                  ENDIF
+
+                  ! DIC straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5))
+
+                  ! Alkalinity straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3a) = outincs(ji,jj,jk,i_tracer(6))
+      
+               END DO
+            END DO
+         END DO
+      
+      ENDIF
       !
-      ! Now set MLD to bottom of a level and propagate incs equally through mixed layer
+      ! Now set MLD to bottom of a level and propagate chlorophyll incs equally through mixed layer
+      ! If balancing, should really relate this back to phytoplankton, but stick with this for now
       !
       DO jj = 1, jpj
@@ -201,6 +549,4 @@
       END DO
 
-      ! Multivariate balancing forthcoming...
-
    END SUBROUTINE asm_logchl_bal_ersem
 
