Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_ersem.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90	(revision 10574)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_ersem.F90	(revision 10622)
@@ -1,6 +1,6 @@
-MODULE asmphyto2dbal_medusa
+MODULE asmphyto2dbal_ersem
    !!======================================================================
-   !!                       ***  MODULE asmphyto2dbal_medusa  ***
-   !! Calculate increments to MEDUSA based on surface phyto2d increments
+   !!                       ***  MODULE asmphyto2dbal_ersem  ***
+   !! Calculate increments to ERSEM based on surface phyto2d increments
    !!
    !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
@@ -10,21 +10,18 @@
    !!
    !!======================================================================
-   !! History :  3.6  ! 2017-08 (D. Ford)  Adapted from asmphyto2dbal_hadocc
+   !! History :  3.6  ! 2019-01 (D. Ford)  Adapted from asmphyto2dbal_medusa
    !!----------------------------------------------------------------------
-#if defined key_asminc && defined key_medusa
+#if defined key_asminc && defined key_fabm
    !!----------------------------------------------------------------------
    !! 'key_asminc'          : assimilation increment interface
-   !! 'key_medusa'          : MEDUSA model
+   !! 'key_fabm'            : FABM-ERSEM model
    !!----------------------------------------------------------------------
-   !! asm_phyto2d_bal_medusa : routine to calculate increments to MEDUSA
+   !! asm_phyto2d_bal_ersem : routine to calculate increments to ERSEM
    !!----------------------------------------------------------------------
    USE par_kind,      ONLY: wp             ! kind parameters
    USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
    USE dom_oce,       ONLY: gdepw_n        ! domain information
-   USE zdftmx,        ONLY: ln_tmx_itf, &  ! Indonesian Throughflow
-      &                     mask_itf       ! tidal mixing mask
    USE iom                                 ! i/o
-   USE sms_medusa                          ! MEDUSA parameters
-   USE par_medusa                          ! MEDUSA parameters
+   USE par_fabm                            ! FABM-ERSEM parameters
    USE par_trc,       ONLY: jptra          ! Tracer parameters
    USE bioanalysis                         ! Nitrogen balancing
@@ -33,5 +30,5 @@
    PRIVATE                   
 
-   PUBLIC asm_phyto2d_bal_medusa
+   PUBLIC asm_phyto2d_bal_ersem
 
    ! Default values for biological assimilation parameters
@@ -68,29 +65,31 @@
 CONTAINS
 
-   SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      &
+   SUBROUTINE asm_phyto2d_bal_ersem(  ld_chltot,                      &
       &                               pinc_chltot,                    &
       &                               ld_chldia,                      &
       &                               pinc_chldia,                    &
-      &                               ld_chlnon,                      &
-      &                               pinc_chlnon,                    &
-      &                               ld_phytot,                      &
-      &                               pinc_phytot,                    &
-      &                               ld_phydia,                      &
-      &                               pinc_phydia,                    &
-      &                               ld_phynon,                      &
-      &                               pinc_phynon,                    &
+      &                               ld_chlnan,                      &
+      &                               pinc_chlnan,                    &
+      &                               ld_chlpic,                      &
+      &                               pinc_chlpic,                    &
+      &                               ld_chldin,                      &
+      &                               pinc_chldin,                    &
       &                               pincper,                        &
       &                               p_maxchlinc, ld_phytobal, pmld, &
       &                               pgrow_avg_bkg, ploss_avg_bkg,   &
       &                               phyt_avg_bkg, mld_max_bkg,      &
+      &                               totalk_bkg,                     &
       &                               tracer_bkg, phyto2d_balinc )
       !!---------------------------------------------------------------------------
-      !!                    ***  ROUTINE asm_phyto2d_bal_medusa  ***
+      !!                    ***  ROUTINE asm_phyto2d_bal_ersem  ***
       !!
-      !! ** Purpose :   calculate increments to MEDUSA from 2d phytoplankton increments
+      !! ** Purpose :   calculate increments to ERSEM from 2d phytoplankton increments
       !!
-      !! ** Method  :   average up MEDUSA to look like HadOCC
+      !! ** Method  :   EITHER (ld_phytobal == .TRUE.):
+      !!                average up ERSEM to look like HadOCC
       !!                call nitrogen balancing scheme
       !!                separate back out to MEDUSA
+      !!                OR (ld_phytobal == .FALSE.):
+      !!                calculate increments to maintain background stoichiometry
       !!
       !! ** Action  :   populate phyto2d_balinc
@@ -98,4 +97,5 @@
       !! References :   Hemmings et al., 2008, J. Mar. Res.
       !!                Ford et al., 2012, Ocean Sci.
+      !!                Skakala et al., 2018, JGR
       !!---------------------------------------------------------------------------
       !!
@@ -104,12 +104,10 @@
       LOGICAL,  INTENT(in   )                               :: ld_chldia      ! Assim chldia y/n
       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldia    ! chldia increments
-      LOGICAL,  INTENT(in   )                               :: ld_chlnon      ! Assim chlnon y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlnon    ! chlnon increments
-      LOGICAL,  INTENT(in   )                               :: ld_phytot      ! Assim phytot y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phytot    ! phytot increments
-      LOGICAL,  INTENT(in   )                               :: ld_phydia      ! Assim phydia y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phydia    ! phydia increments
-      LOGICAL,  INTENT(in   )                               :: ld_phynon      ! Assim phynon y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phynon    ! phynon increments
+      LOGICAL,  INTENT(in   )                               :: ld_chlnan      ! Assim chlnan y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlnan    ! chlnan increments
+      LOGICAL,  INTENT(in   )                               :: ld_chlpic      ! Assim chlpic y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlpic    ! chlpic increments
+      LOGICAL,  INTENT(in   )                               :: ld_chldin      ! Assim chldin y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldin    ! chldin increments
       REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
       REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
@@ -120,4 +118,5 @@
       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)       :: totalk_bkg     ! Total alkalinity
       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
@@ -126,21 +125,66 @@
       INTEGER                                               :: jkmax          ! Loop index
       INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
+      REAL(wp)                                              :: zmassc         ! Carbon molar mass
+      REAL(wp)                                              :: zmassn         ! Nitrogen molar mass
+      REAL(wp)                                              :: z4qnc          ! Z4/qnc (mesozoo N:C)
       REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
       REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
       REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
-      REAL(wp)                                              :: zfrac          ! Fraction
-      REAL(wp)                                              :: zfrac_chn      ! Fraction of jpchn
-      REAL(wp)                                              :: zfrac_chd      ! Fraction of jpchd
-      REAL(wp)                                              :: zfrac_phn      ! Fraction of jpphn
-      REAL(wp)                                              :: zfrac_phd      ! Fraction of jpphd
-      REAL(wp)                                              :: zfrac_zmi      ! Fraction of jpzmi
-      REAL(wp)                                              :: zfrac_zme      ! Fraction of jpzme
-      REAL(wp)                                              :: zrat_pds_phd   ! Ratio of jppds:jpphd
-      REAL(wp)                                              :: zrat_chd_phd   ! Ratio of jpchd:jpphd
-      REAL(wp)                                              :: zrat_chn_phn   ! Ratio of jpchn:jpphn
-      REAL(wp)                                              :: zrat_phn_chn   ! Ratio of jpphn:jpchn
-      REAL(wp)                                              :: zrat_phd_chd   ! Ratio of jpphd:jpchd
-      REAL(wp)                                              :: zrat_pds_chd   ! Ratio of jppds:jpchd
-      REAL(wp)                                              :: zrat_dtc_det   ! Ratio of jpdtc:jpdet
+      REAL(wp)                                              :: zfrac          ! Fractions
+      REAL(wp)                                              :: zfrac_chl1     ! 
+      REAL(wp)                                              :: zfrac_chl2     ! 
+      REAL(wp)                                              :: zfrac_chl3     ! 
+      REAL(wp)                                              :: zfrac_chl4     ! 
+      REAL(wp)                                              :: zfrac_p1n      ! 
+      REAL(wp)                                              :: zfrac_p2n      ! 
+      REAL(wp)                                              :: zfrac_p3n      ! 
+      REAL(wp)                                              :: zfrac_p4n      ! 
+      REAL(wp)                                              :: zfrac_z4n      ! 
+      REAL(wp)                                              :: zfrac_z5n      ! 
+      REAL(wp)                                              :: zfrac_z6n      ! 
+      REAL(wp)                                              :: zfrac_n3n      ! 
+      REAL(wp)                                              :: zfrac_n4n      ! 
+      REAL(wp)                                              :: zfrac_r4n      ! 
+      REAL(wp)                                              :: zfrac_r6n      ! 
+      REAL(wp)                                              :: zfrac_r8n      ! 
+      REAL(wp)                                              :: zrat_chl1_p1n  ! Ratios
+      REAL(wp)                                              :: zrat_p1c_p1n   ! 
+      REAL(wp)                                              :: zrat_p1p_p1n   ! 
+      REAL(wp)                                              :: zrat_p1s_p1n   ! 
+      REAL(wp)                                              :: zrat_chl2_p2n  ! 
+      REAL(wp)                                              :: zrat_p2c_p2n   ! 
+      REAL(wp)                                              :: zrat_p2p_p2n   ! 
+      REAL(wp)                                              :: zrat_chl3_p3n  ! 
+      REAL(wp)                                              :: zrat_p3c_p3n   ! 
+      REAL(wp)                                              :: zrat_p3p_p3n   ! 
+      REAL(wp)                                              :: zrat_chl4_p4n  ! 
+      REAL(wp)                                              :: zrat_p4c_p4n   ! 
+      REAL(wp)                                              :: zrat_p4p_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_r4c_r4n   ! 
+      REAL(wp)                                              :: zrat_r4p_r4n   ! 
+      REAL(wp)                                              :: zrat_r6c_r6n   ! 
+      REAL(wp)                                              :: zrat_r6p_r6n   ! 
+      REAL(wp)                                              :: zrat_r6s_r6n   ! 
+      REAL(wp)                                              :: zrat_r8c_r8n   ! 
+      REAL(wp)                                              :: zrat_r8p_r8n   ! 
+      REAL(wp)                                              :: zrat_r8s_r8n   ! 
+      REAL(wp)                                              :: zrat_p1c_chl1  ! 
+      REAL(wp)                                              :: zrat_p1n_chl1  ! 
+      REAL(wp)                                              :: zrat_p1p_chl1  ! 
+      REAL(wp)                                              :: zrat_p1s_chl1  ! 
+      REAL(wp)                                              :: zrat_p2c_chl2  ! 
+      REAL(wp)                                              :: zrat_p2n_chl2  ! 
+      REAL(wp)                                              :: zrat_p2p_chl2  ! 
+      REAL(wp)                                              :: zrat_p3c_chl3  ! 
+      REAL(wp)                                              :: zrat_p3n_chl3  ! 
+      REAL(wp)                                              :: zrat_p3p_chl3  ! 
+      REAL(wp)                                              :: zrat_p4c_chl4  ! 
+      REAL(wp)                                              :: zrat_p4n_chl4  ! 
+      REAL(wp)                                              :: zrat_p4p_chl4  ! 
       REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
       REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
@@ -150,5 +194,14 @@
       REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
       REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_temp
       !!---------------------------------------------------------------------------
+      
+      ! Set parameters
+      zmassc = 12.01
+      zmassn = 14.01
+      z4qnc  = 0.0126
+      !z4qnc  = model%state_variables(jp_fabm_z4c)%parameters%qnc%value
+      !z4qnc  = get_property_by_name(model%state_variables(jp_fabm_z4c)%parameters, 'qnc')
+      IF (lwp) WRITE(numout,*) 'z4qnc = ', z4qnc
 
       ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
@@ -160,36 +213,66 @@
                END DO
             END DO
-         ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
+         ELSE
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  pinc_chltot(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnon(ji,jj)
-                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
-                  IF ( pinc_chltot(ji,jj) .NE. ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) ) THEN
-                     zfrac = pinc_chltot(ji,jj) / ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) )
-                     pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac
-                     pinc_chlnon(ji,jj) = pinc_chlnon(ji,jj) * zfrac
-                  ENDIF
-               END DO
-            END DO
-         ELSE IF ( ld_chldia ) THEN
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pinc_chldia(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chldia(ji,jj), p_maxchlinc ) )
-                  pinc_chltot(ji,jj) = pinc_chldia(ji,jj)
-               END DO
-            END DO
-         ELSE IF ( ld_chlnon ) THEN
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pinc_chlnon(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chlnon(ji,jj), p_maxchlinc ) )
-                  pinc_chltot(ji,jj) = pinc_chlnon(ji,jj)
+                  IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
+                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
+                        &                      pinc_chlpic(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
+                        &                      pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlpic .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + &
+                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + &
+                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlnan ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chlpic ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlpic(ji,jj)
+                  ELSE IF ( ld_chldia .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chlnan .AND. ld_chlpic ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chlpic(ji,jj)
+                  ELSE IF ( ld_chlnan .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chlpic .AND. ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
+                  ELSE IF ( ld_chldia ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj)
+                  ELSE IF ( ld_chlnan ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj)
+                  ELSE IF ( ld_chlpic ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj)
+                  ELSE IF ( ld_chldin ) THEN
+                     pinc_chltot_temp(ji,jj) = pinc_chldin(ji,jj)
+                  ENDIF
+                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_temp(ji,jj), p_maxchlinc ) )
+                  IF ( pinc_chltot(ji,jj) .NE. pinc_chltot_temp(ji,jj) ) THEN
+                     zfrac = pinc_chltot(ji,jj) / pinc_chltot_temp(ji,jj)
+                     IF ( ld_chldia ) THEN
+                        pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac
+                     ENDIF
+                     IF ( ld_chlnan ) THEN
+                        pinc_chlnan(ji,jj) = pinc_chlnan(ji,jj) * zfrac
+                     ENDIF
+                     IF ( ld_chlpic ) THEN
+                        pinc_chlpic(ji,jj) = pinc_chlpic(ji,jj) * zfrac
+                     ENDIF
+                     IF ( ld_chldin ) THEN
+                        pinc_chldin(ji,jj) = pinc_chldin(ji,jj) * zfrac
+                     ENDIF
+                  ENDIF
                END DO
             END DO
          ENDIF
       ENDIF
-
-      IF ( ld_phytot .OR. ld_phydia .OR. ld_phynon ) THEN
-         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
-      ENDIF
+      
+      ! Initialise balancing increments
+      phyto2d_balinc(:,:,:,:) = 0.0
       
       IF ( ld_phytobal ) THEN   ! Nitrogen balancing
@@ -197,6 +280,6 @@
          ! 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
+         ! 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.
@@ -211,7 +294,7 @@
          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(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
@@ -250,27 +333,40 @@
 
          ! 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)
+         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_fabm_n3n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_n4n)
+         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_fabm_p1n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_p2n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_p3n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_p4n)
+         bstate(:,:,:,i_tracer(3)) = (tracer_bkg(:,:,:,jp_fabm_z4c) * z4qnc) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_z5n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_z6n)
+         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_fabm_r4n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_r6n) + &
+            &                        tracer_bkg(:,:,:,jp_fabm_r8n)
+         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_fabm_o3c)
+         bstate(:,:,:,i_tracer(6)) = totalk_bkg(:,:,:)
 
          ! 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
+         ! and nitrogen to biomass equivalent for PZD (hardwire as per HadOCC)
          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 ) )
+               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                  &   tracer_bkg(ji,jj,1,jp_fabm_chl3) + tracer_bkg(ji,jj,1,jp_fabm_chl4) ) .GT. 0.0 ) THEN
+                  cchl_p(ji,jj) = zmassc * ( tracer_bkg(ji,jj,1,jp_fabm_p1c) +      &
+                     &                       tracer_bkg(ji,jj,1,jp_fabm_p2c) +      &
+                     &                       tracer_bkg(ji,jj,1,jp_fabm_p3c) +      &
+                     &                       tracer_bkg(ji,jj,1,jp_fabm_p4c)   ) /  &
+                     &            ( tracer_bkg(ji,jj,1,jp_fabm_chl1) +             &
+                     &              tracer_bkg(ji,jj,1,jp_fabm_chl2) +             &
+                     &              tracer_bkg(ji,jj,1,jp_fabm_chl3) +             &
+                     &              tracer_bkg(ji,jj,1,jp_fabm_chl4)   )
                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 )
+         n2be_p = zmassn + ( zmassc * modparm(10) )
+         n2be_z = zmassn + ( zmassc * modparm(11) )
+         n2be_d = zmassn + ( zmassc * modparm(12) )
 
          ! Call nitrogen balancing routine
@@ -288,5 +384,4 @@
          
          ! Loop over each grid point partioning the increments
-         phyto2d_balinc(:,:,:,:) = 0.0
          DO jk = 1, jpk
             DO jj = 1, jpj
@@ -294,71 +389,180 @@
 
                   ! Phytoplankton
-                  IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. &
-                     & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. &
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p2n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p3n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p4n) > 0.0 ) .AND. &
                      & ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN
                      IF ( ld_chltot ) THEN
                         ! Phytoplankton nitrogen split up based on existing ratios
-                        zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / &
-                           &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
-                        zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / &
-                           &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
-                     ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
+                        zfrac_p1n = tracer_bkg(ji,jj,jk,jp_fabm_p1n) / &
+                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) )
+                        zfrac_p2n = tracer_bkg(ji,jj,jk,jp_fabm_p2n) / &
+                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) )
+                        zfrac_p3n = tracer_bkg(ji,jj,jk,jp_fabm_p3n) / &
+                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) )
+                        zfrac_p4n = tracer_bkg(ji,jj,jk,jp_fabm_p4n) / &
+                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + &
+                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) )
+                     ELSE
                         ! Phytoplankton nitrogen split up based on assimilation increments
-                        zfrac_phn = pinc_chlnon(ji,jj) / pinc_chltot(ji,jj)
-                        zfrac_phd = pinc_chldia(ji,jj) / pinc_chltot(ji,jj)
+                        zfrac_p1n = pinc_chldia(ji,jj) / pinc_chltot(ji,jj)
+                        zfrac_p2n = pinc_chlnan(ji,jj) / pinc_chltot(ji,jj)
+                        zfrac_p3n = pinc_chlpic(ji,jj) / pinc_chltot(ji,jj)
+                        zfrac_p4n = pinc_chldin(ji,jj) / pinc_chltot(ji,jj)
                      ENDIF
-
-                     ! Phytoplankton silicate split up based on existing ratios
-                     zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd)
                      
-                     ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen
-                     ! Not using pinc_chltot directly as it's only 2D
-                     ! This method should give same results at surface as splitting pinc_chltot would
-                     zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn)
-                     zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd)
+                     ! Other phytoplankton variables split up based on existing ratios with nitrogen
+                     zrat_chl1_p1n = tracer_bkg(ji,jj,jk,jp_fabm_chl1) / tracer_bkg(ji,jj,jk,jp_fabm_p1n)
+                     zrat_p1c_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1c)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n)
+                     zrat_p1p_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1p)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n)
+                     zrat_p1s_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1s)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n)
+                     zrat_chl2_p2n = tracer_bkg(ji,jj,jk,jp_fabm_chl2) / tracer_bkg(ji,jj,jk,jp_fabm_p2n)
+                     zrat_p2c_p2n  = tracer_bkg(ji,jj,jk,jp_fabm_p2c)  / tracer_bkg(ji,jj,jk,jp_fabm_p2n)
+                     zrat_p2p_p2n  = tracer_bkg(ji,jj,jk,jp_fabm_p2p)  / tracer_bkg(ji,jj,jk,jp_fabm_p2n)
+                     zrat_chl3_p3n = tracer_bkg(ji,jj,jk,jp_fabm_chl3) / tracer_bkg(ji,jj,jk,jp_fabm_p3n)
+                     zrat_p3c_p3n  = tracer_bkg(ji,jj,jk,jp_fabm_p3c)  / tracer_bkg(ji,jj,jk,jp_fabm_p3n)
+                     zrat_p3p_p3n  = tracer_bkg(ji,jj,jk,jp_fabm_p3p)  / tracer_bkg(ji,jj,jk,jp_fabm_p3n)
+                     zrat_chl4_p4n = tracer_bkg(ji,jj,jk,jp_fabm_chl4) / tracer_bkg(ji,jj,jk,jp_fabm_p4n)
+                     zrat_p4c_p4n  = tracer_bkg(ji,jj,jk,jp_fabm_p4c)  / tracer_bkg(ji,jj,jk,jp_fabm_p4n)
+                     zrat_p4p_p4n  = tracer_bkg(ji,jj,jk,jp_fabm_p4p)  / tracer_bkg(ji,jj,jk,jp_fabm_p4n)
                      
-                     phyto2d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn
-                     phyto2d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd
-                     phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
-                     phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
-                     phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_chd_phd
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl1) = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_chl1_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1c_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1p_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1s)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1s_p1n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl2) = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_chl2_p2n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_p2c_p2n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_p2p_p2n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl3) = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_chl3_p3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_p3c_p3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_p3p_p3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl4) = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_chl4_p4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_p4c_p4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_p4p_p4n
                   ENDIF
 
                   ! Zooplankton nitrogen split up based on existing ratios
-                  IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN
-                     zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / &
-                        &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
-                     zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / &
-                        &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
-                     phyto2d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi
-                     phyto2d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme
-                  ENDIF
-
-                  ! Nitrogen nutrient straight from balancing scheme
-                  phyto2d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1))
-
-                  ! Nitrogen detritus straight from balancing scheme
-                  phyto2d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4))
+                  ! Update carbon and phosphorus according to existing ratios
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_z5n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_z6n) > 0.0 ) ) THEN
+                     zfrac_z4n = ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) / &
+                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) )
+                     zfrac_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5n) / &
+                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) )
+                     zfrac_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6n) / &
+                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + &
+                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) )
+                     zrat_z4c_z4n = 1.0 / z4qnc
+                     zrat_z5c_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5c) / tracer_bkg(ji,jj,jk,jp_fabm_z5n)
+                     zrat_z5p_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5p) / tracer_bkg(ji,jj,jk,jp_fabm_z5n)
+                     zrat_z6c_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6c) / tracer_bkg(ji,jj,jk,jp_fabm_z6n)
+                     zrat_z6p_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6p) / tracer_bkg(ji,jj,jk,jp_fabm_z6n)
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5c) = phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) * zrat_z5c_z5n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6c) = phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) * zrat_z6c_z6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5p) = phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) * zrat_z5p_z5n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6p) = phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) * zrat_z6p_z6n
+                  ENDIF
+
+                  ! Nitrogen nutrient split between nitrate and ammonium based on existing ratios
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n3n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_n4n) > 0.0 ) ) THEN
+                     zfrac_n3n = tracer_bkg(ji,jj,jk,jp_fabm_n3n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_n4n))
+                     zfrac_n4n = tracer_bkg(ji,jj,jk,jp_fabm_n4n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_n4n))
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n
+                  ENDIF
+
+                  ! Detritus nitrogen split up based on existing ratios
+                  ! Update carbon and phosphorus according to existing ratios
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_r4n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_r6n) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jp_fabm_r8n) > 0.0 ) ) THEN
+                     zfrac_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n))
+                     zfrac_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n))
+                     zfrac_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8n) / &
+                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + &
+                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n))
+                     zrat_r4c_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4c) / tracer_bkg(ji,jj,jk,jp_fabm_r4n)
+                     zrat_r4p_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4p) / tracer_bkg(ji,jj,jk,jp_fabm_r4n)
+                     zrat_r6c_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6c) / tracer_bkg(ji,jj,jk,jp_fabm_r6n)
+                     zrat_r6p_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6p) / tracer_bkg(ji,jj,jk,jp_fabm_r6n)
+                     zrat_r6s_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6s) / tracer_bkg(ji,jj,jk,jp_fabm_r6n)
+                     zrat_r8c_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8c) / tracer_bkg(ji,jj,jk,jp_fabm_r8n)
+                     zrat_r8p_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8p) / tracer_bkg(ji,jj,jk,jp_fabm_r8n)
+                     zrat_r8s_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8s) / tracer_bkg(ji,jj,jk,jp_fabm_r8n)
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) * zrat_r4c_r4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) * zrat_r4p_r4n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6c_r6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6p_r6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6s) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6s_r6n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8c_r8n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8p_r8n
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8s) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8s_r8n
+                  ENDIF
 
                   ! DIC straight from balancing scheme
-                  phyto2d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5))
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5))
 
                   ! Alkalinity straight from balancing scheme
-                  phyto2d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6))
-
-                  ! Remove diatom silicate increment from nutrient silicate to conserve mass
-                  IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto2d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
-                     phyto2d_balinc(ji,jj,jk,jpsil) = phyto2d_balinc(ji,jj,jk,jppds) * (-1.0)
-                  ENDIF
-
-                  ! Carbon detritus based on existing ratios
-                  IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN
-                     zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet)
-                     phyto2d_balinc(ji,jj,jk,jpdtc) = phyto2d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
-                  ENDIF
-
-                  ! Do nothing with iron or oxygen for the time being
-                  phyto2d_balinc(ji,jj,jk,jpfer) = 0.0
-                  phyto2d_balinc(ji,jj,jk,jpoxy) = 0.0
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_o3ba) = outincs(ji,jj,jk,i_tracer(6))
+
+                  ! Remove P/R silicon increments from silicate to conserve mass
+                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_p1s) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r6s) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r8s)
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n5s) - zfrac ) > 0.0 ) THEN
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_n5s) = zfrac * (-1.0)
+                  ENDIF
+
+                  ! Remove P/Z/R phosphorus increments from phosphate to conserve mass
+                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_p1p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p2p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p3p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p4p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_z5p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_z6p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r4p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r6p) + &
+                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r8p)
+                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n1p) - zfrac ) > 0.0 ) THEN
+                     phyto2d_balinc(ji,jj,jk,jp_fabm_n1p) = zfrac * (-1.0)
+                  ENDIF
                   
                END DO
@@ -366,41 +570,87 @@
          END DO
       
-      ELSE   ! No nitrogen balancing
-      
-         ! Initialise individual chlorophyll increments to zero
-         phyto2d_balinc(:,:,:,jpchn) = 0.0
-         phyto2d_balinc(:,:,:,jpchd) = 0.0
+      ELSE   ! No nitrogen balancing - just update phytoplankton
          
          ! Split up total surface chlorophyll increments
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( ( tracer_bkg(ji,jj,1,jpchn) > 0.0 ) .AND. &
-                  & ( tracer_bkg(ji,jj,1,jpchd) > 0.0 ) ) THEN
+               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_chl1) > 0.0 ) .AND. &
+                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl2) > 0.0 ) .AND. &
+                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl3) > 0.0 ) .AND. &
+                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl4) > 0.0 ) ) THEN
                   IF ( ld_chltot ) THEN
                      ! Chlorophyll split up based on existing ratios
-                     zfrac_chn = tracer_bkg(ji,jj,1,jpchn) / &
-                        &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) )
-                     zfrac_chd = tracer_bkg(ji,jj,1,jpchd) / &
-                        &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) )
-                     phyto2d_balinc(ji,jj,1,jpchn) = pinc_chltot(ji,jj) * zfrac_chn
-                     phyto2d_balinc(ji,jj,1,jpchd) = pinc_chltot(ji,jj) * zfrac_chd
+                     zfrac_chl1 = tracer_bkg(ji,jj,1,jp_fabm_chl1) /  &
+                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) )
+                     zfrac_chl2 = tracer_bkg(ji,jj,1,jp_fabm_chl2) /  &
+                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) )
+                     zfrac_chl3 = tracer_bkg(ji,jj,1,jp_fabm_chl3) /  &
+                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) )
+                     zfrac_chl4 = tracer_bkg(ji,jj,1,jp_fabm_chl4) /  &
+                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + &
+                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) )
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl1) = pinc_chltot(ji,jj) * zfrac_chl1
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl2) = pinc_chltot(ji,jj) * zfrac_chl2
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl3) = pinc_chltot(ji,jj) * zfrac_chl3
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl4) = pinc_chltot(ji,jj) * zfrac_chl4
                   ENDIF
                   IF( ld_chldia ) THEN
-                     phyto2d_balinc(ji,jj,1,jpchd) = pinc_chldia(ji,jj)
-                  ENDIF
-                  IF( ld_chlnon ) THEN
-                     phyto2d_balinc(ji,jj,1,jpchn) = pinc_chlnon(ji,jj)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl1) = pinc_chldia(ji,jj)
+                  ENDIF
+                  IF( ld_chlnan ) THEN
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl2) = pinc_chlnan(ji,jj)
+                  ENDIF
+                  IF( ld_chlpic ) THEN
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl3) = pinc_chlpic(ji,jj)
+                  ENDIF
+                  IF( ld_chldin ) THEN
+                     phyto2d_balinc(ji,jj,1,jp_fabm_chl4) = pinc_chldin(ji,jj)
                   ENDIF
                   
-                  ! Maintain stoichiometric ratios of nitrogen and silicate
-                  IF ( ld_chltot .OR. ld_chlnon ) THEN
-                     zrat_phn_chn = tracer_bkg(ji,jj,1,jpphn) / tracer_bkg(ji,jj,1,jpchn)
-                     phyto2d_balinc(ji,jj,1,jpphn) = phyto2d_balinc(ji,jj,1,jpchn) * zrat_phn_chn
-                  ENDIF
+                  ! Maintain stoichiometric ratios of carbon, nitrogen, phosphorus and silicon
                   IF ( ld_chltot .OR. ld_chldia ) THEN
-                     zrat_phd_chd = tracer_bkg(ji,jj,1,jpphd) / tracer_bkg(ji,jj,1,jpchd)
-                     phyto2d_balinc(ji,jj,1,jpphd) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_phd_chd
-                     zrat_pds_chd = tracer_bkg(ji,jj,1,jppds) / tracer_bkg(ji,jj,1,jpchd)
-                     phyto2d_balinc(ji,jj,1,jppds) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_pds_chd
+                     zrat_p1c_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1c) / tracer_bkg(ji,jj,1,jp_fabm_chl1)
+                     zrat_p1n_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1n) / tracer_bkg(ji,jj,1,jp_fabm_chl1)
+                     zrat_p1p_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1p) / tracer_bkg(ji,jj,1,jp_fabm_chl1)
+                     zrat_p1s_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1s) / tracer_bkg(ji,jj,1,jp_fabm_chl1)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p1c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1c_chl1
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p1n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1n_chl1
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p1p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1p_chl1
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p1s) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1s_chl1
+                  ENDIF
+                  IF ( ld_chltot .OR. ld_chlnan ) THEN
+                     zrat_p2c_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2c) / tracer_bkg(ji,jj,1,jp_fabm_chl2)
+                     zrat_p2n_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2n) / tracer_bkg(ji,jj,1,jp_fabm_chl2)
+                     zrat_p2p_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2p) / tracer_bkg(ji,jj,1,jp_fabm_chl2)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p2c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2c_chl2
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p2n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2n_chl2
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p2p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2p_chl2
+                  ENDIF
+                  IF ( ld_chltot .OR. ld_chlpic ) THEN
+                     zrat_p3c_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3c) / tracer_bkg(ji,jj,1,jp_fabm_chl3)
+                     zrat_p3n_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3n) / tracer_bkg(ji,jj,1,jp_fabm_chl3)
+                     zrat_p3p_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3p) / tracer_bkg(ji,jj,1,jp_fabm_chl3)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p3c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3c_chl3
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p3n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3n_chl3
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p3p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3p_chl3
+                  ENDIF
+                  IF ( ld_chltot .OR. ld_chldin ) THEN
+                     zrat_p4c_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4c) / tracer_bkg(ji,jj,1,jp_fabm_chl4)
+                     zrat_p4n_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4n) / tracer_bkg(ji,jj,1,jp_fabm_chl4)
+                     zrat_p4p_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4p) / tracer_bkg(ji,jj,1,jp_fabm_chl4)
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p4c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4c_chl4
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p4n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4n_chl4
+                     phyto2d_balinc(ji,jj,1,jp_fabm_p4p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4p_chl4
                   ENDIF
                ENDIF
@@ -422,9 +672,21 @@
                !
                DO jk = 2, jkmax
-                  phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,1,jpchn)
-                  phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,1,jpchd)
-                  phyto2d_balinc(ji,jj,jk,jpphn) = phyto2d_balinc(ji,jj,1,jpphn)
-                  phyto2d_balinc(ji,jj,jk,jpphd) = phyto2d_balinc(ji,jj,1,jpphd)
-                  phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,1,jppds)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl1) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1c)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1n)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1p)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1s)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1s)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl2) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2c)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2n)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2p)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl3) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3c)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3n)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3p)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl4) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4c)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4n)
+                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4p)
                END DO
                !
@@ -432,30 +694,7 @@
          END DO
 
-         ! Set other balancing increments to zero
-         phyto2d_balinc(:,:,:,jpzmi) = 0.0
-         phyto2d_balinc(:,:,:,jpzme) = 0.0
-         phyto2d_balinc(:,:,:,jpdin) = 0.0
-         phyto2d_balinc(:,:,:,jpsil) = 0.0
-         phyto2d_balinc(:,:,:,jpfer) = 0.0
-         phyto2d_balinc(:,:,:,jpdet) = 0.0
-         phyto2d_balinc(:,:,:,jpdtc) = 0.0
-         phyto2d_balinc(:,:,:,jpdic) = 0.0
-         phyto2d_balinc(:,:,:,jpalk) = 0.0
-         phyto2d_balinc(:,:,:,jpoxy) = 0.0
-
       ENDIF
-      
-      ! If performing extra tidal mixing in the Indonesian Throughflow,
-      ! increments have been found to make the carbon cycle unstable
-      ! Therefore, mask these out
-      IF ( ln_tmx_itf ) THEN
-         DO jn = 1, jptra
-            DO jk = 1, jpk
-               phyto2d_balinc(:,:,jk,jn) = phyto2d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
-            END DO
-         END DO
-      ENDIF
-
-   END SUBROUTINE asm_phyto2d_bal_medusa
+
+   END SUBROUTINE asm_phyto2d_bal_ersem
 
 #else
@@ -464,33 +703,30 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      &
-      &                              pinc_chltot,                    &
-      &                              ld_chldia,                      &
-      &                              pinc_chldia,                    &
-      &                              ld_chlnon,                      &
-      &                              pinc_chlnon,                    &
-      &                              ld_phytot,                      &
-      &                              pinc_phytot,                    &
-      &                              ld_phydia,                      &
-      &                              pinc_phydia,                    &
-      &                              ld_phynon,                      &
-      &                              pinc_phynon,                    &
-      &                              pincper,                        &
-      &                              p_maxchlinc, ld_phytobal, pmld, &
-      &                              pgrow_avg_bkg, ploss_avg_bkg,   &
-      &                              phyt_avg_bkg, mld_max_bkg,      &
-      &                              tracer_bkg, phyto2d_balinc )
+   SUBROUTINE asm_phyto2d_bal_ersem(  ld_chltot,                      &
+      &                               pinc_chltot,                    &
+      &                               ld_chldia,                      &
+      &                               pinc_chldia,                    &
+      &                               ld_chlnan,                      &
+      &                               pinc_chlnan,                    &
+      &                               ld_chlpic,                      &
+      &                               pinc_chlpic,                    &
+      &                               ld_chldin,                      &
+      &                               pinc_chldin,                    &
+      &                               pincper,                        &
+      &                               p_maxchlinc, ld_phytobal, pmld, &
+      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                               phyt_avg_bkg, mld_max_bkg,      &
+      &                               totalk_bkg,                     &
+      &                               tracer_bkg, phyto2d_balinc )
       LOGICAL :: ld_chltot
       REAL    :: pinc_chltot(:,:)
       LOGICAL :: ld_chldia
       REAL    :: pinc_chldia(:,:)
-      LOGICAL :: ld_chlnon
-      REAL    :: pinc_chlnon(:,:)
-      LOGICAL :: ld_phytot
-      REAL    :: pinc_phytot(:,:)
-      LOGICAL :: ld_phydia
-      REAL    :: pinc_phydia(:,:)
-      LOGICAL :: ld_phynon
-      REAL    :: pinc_phynon(:,:)
+      LOGICAL :: ld_chlnan
+      REAL    :: pinc_chlnan(:,:)
+      LOGICAL :: ld_chlpic
+      REAL    :: pinc_chlpic(:,:)
+      LOGICAL :: ld_chldin
+      REAL    :: pinc_chldin(:,:)
       REAL    :: pincper
       REAL    :: p_maxchlinc
@@ -501,10 +737,11 @@
       REAL    :: phyt_avg_bkg(:,:)
       REAL    :: mld_max_bkg(:,:)
+      REAL    :: totalk_bkg(:,:,:)
       REAL    :: tracer_bkg(:,:,:,:)
       REAL    :: phyto2d_balinc(:,:,:,:)
-      WRITE(*,*) 'asm_phyto2d_bal_medusa: You should not have seen this print! error?'
-   END SUBROUTINE asm_phyto2d_bal_medusa
+      WRITE(*,*) 'asm_phyto2d_bal_ersem: You should not have seen this print! error?'
+   END SUBROUTINE asm_phyto2d_bal_ersem
 #endif
 
    !!======================================================================
-END MODULE asmphyto2dbal_medusa
+END MODULE asmphyto2dbal_ersem
