Index: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
===================================================================
--- branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 13191)
+++ branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 13354)
@@ -18,5 +18,7 @@
    !! asm_bgc_bal_wri       : write out bgc balancing increments
    !! asm_bgc_bkg_wri       : write out bgc background
-   !! phyto2d_asm_inc       : apply the ocean colour increments
+   !! asm_bgc_unlog_2d      : calculate non-log versions of 2D log increments
+   !! asm_bgc_unlog_3d      : calculate non-log versions of 3D log increments
+   !! phyto2d_asm_inc       : apply the 2D phytoplankton increments
    !! phyto3d_asm_inc       : apply the 3D phytoplankton increments
    !! pco2_asm_inc          : apply the pCO2/fCO2 increments
@@ -56,6 +58,6 @@
 #endif
 #if defined key_medusa
-   USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
-      & asm_phyto2d_bal_medusa
+   USE asmphytobal_medusa, ONLY: & ! phytoplankton balancing for MEDUSA
+      & asm_phyto_bal_medusa
    USE asmpco2bal, ONLY:    & ! pCO2 balancing for MEDUSA
       & asm_pco2_bal
@@ -70,8 +72,11 @@
       & ploss_avg,          &
       & phyt_avg,           &
+      & pgrow_avg_3d,       &
+      & ploss_avg_3d,       &
+      & phyt_avg_3d,        &
       & mld_max
 #elif defined key_hadocc
-   USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC
-      & asm_phyto2d_bal_hadocc
+   USE asmphytobal_hadocc, ONLY: & ! phytoplankton balancing for HadOCC
+      & asm_phyto_bal_hadocc
    USE asmpco2bal, ONLY:    & ! pCO2 balancing for HadOCC
       & asm_pco2_bal
@@ -82,4 +87,7 @@
       & ploss_avg,          &
       & phyt_avg,           &
+      & pgrow_avg_3d,       &
+      & ploss_avg_3d,       &
+      & phyt_avg_3d,        &
       & mld_max,            &
       & HADOCC_CHL
@@ -98,4 +106,6 @@
    PUBLIC  asm_bgc_bal_wri        ! called by nemo_gcm in nemogcm.F90
    PUBLIC  asm_bgc_bkg_wri        ! called by asm_bkg_wri in asmbkg.F90
+   PRIVATE asm_bgc_unlog_2d       ! called by phyto2d_asm_inc
+   PRIVATE asm_bgc_unlog_3d       ! called by phyto3d_asm_inc
    PUBLIC  phyto2d_asm_inc        ! called by bgc_asm_inc in asminc.F90
    PUBLIC  phyto3d_asm_inc        ! called by bgc_asm_inc in asminc.F90
@@ -168,4 +178,7 @@
    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg   ! Background phyto loss
    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg    ! Background phyto conc
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pgrow_avg_3d_bkg   ! Background phyto growth
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: ploss_avg_3d_bkg   ! Background phyto loss
+   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: phyt_avg_3d_bkg    ! Background phyto conc
    REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg     ! Background max MLD
    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg      ! Background tracer state
@@ -212,7 +225,8 @@
          & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. &
          & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
-         & ( .NOT. ln_slphynoninc ) ) THEN
+         & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
+         & ( .NOT. ln_pchltotinc  ) ) THEN
          CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
-            &           ' if not assimilating ocean colour,',                   &
+            &           ' if not assimilating phytoplankton,',                   &
             &           ' so ln_phytobal will be set to .false.')
          ln_phytobal = .FALSE.
@@ -524,4 +538,7 @@
          ALLOCATE( ploss_avg_bkg(jpi,jpj)        )
          ALLOCATE( phyt_avg_bkg(jpi,jpj)         )
+         ALLOCATE( pgrow_avg_3d_bkg(jpi,jpj,jpk) )
+         ALLOCATE( ploss_avg_3d_bkg(jpi,jpj,jpk) )
+         ALLOCATE( phyt_avg_3d_bkg(jpi,jpj,jpk)  )
          ALLOCATE( mld_max_bkg(jpi,jpj)          )
          ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
@@ -529,4 +546,7 @@
          ploss_avg_bkg(:,:)  = 0.0
          phyt_avg_bkg(:,:)   = 0.0
+         pgrow_avg_3d_bkg(:,:,:)  = 0.0
+         ploss_avg_3d_bkg(:,:,:)  = 0.0
+         phyt_avg_3d_bkg(:,:,:)   = 0.0
          mld_max_bkg(:,:)    = 0.0
          tracer_bkg(:,:,:,:) = 0.0
@@ -564,8 +584,14 @@
             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, 'pgrow_avg_3d', pgrow_avg_3d_bkg )
+            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg_3d', ploss_avg_3d_bkg )
+            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg_3d',  phyt_avg_3d_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)
+            pgrow_avg_3d_bkg(:,:,:) = pgrow_avg_3d_bkg(:,:,:) * tmask(:,:,:)
+            ploss_avg_3d_bkg(:,:,:) = ploss_avg_3d_bkg(:,:,:) * tmask(:,:,:)
+            phyt_avg_3d_bkg(:,:,:)  = phyt_avg_3d_bkg(:,:,:)  * tmask(:,:,:)
             mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1)
 
@@ -726,6 +752,25 @@
             CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) )
             CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) )
+            IF ( ln_phytobal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'phy3d_zmi', phyto3d_balinc(:,:,:,jpzmi) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_zme', phyto3d_balinc(:,:,:,jpzme) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_din', phyto3d_balinc(:,:,:,jpdin) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_sil', phyto3d_balinc(:,:,:,jpsil) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_fer', phyto3d_balinc(:,:,:,jpfer) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_det', phyto3d_balinc(:,:,:,jpdet) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_dtc', phyto3d_balinc(:,:,:,jpdtc) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_dic', phyto3d_balinc(:,:,:,jpdic) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_alk', phyto3d_balinc(:,:,:,jpalk) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_oxy', phyto3d_balinc(:,:,:,jpoxy) )
+            ENDIF
 #elif defined key_hadocc
             CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
+            IF ( ln_phytobal ) THEN
+               CALL iom_rstput( kt, kt, inum, 'phy3d_nut', phyto3d_balinc(:,:,:,jp_had_nut) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_zoo', phyto3d_balinc(:,:,:,jp_had_zoo) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_det', phyto3d_balinc(:,:,:,jp_had_pdn) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_dic', phyto3d_balinc(:,:,:,jp_had_dic) )
+               CALL iom_rstput( kt, kt, inum, 'phy3d_alk', phyto3d_balinc(:,:,:,jp_had_alk) )
+            ENDIF
 #endif
          ENDIF
@@ -792,9 +837,13 @@
       !!------------------------------------------------------------------------
 
-#if defined key_hadocc
+#if defined key_hadocc || defined key_medusa
       CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg             )
       CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg             )
       CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg              )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg_3d', pgrow_avg_3d          )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg_3d', ploss_avg_3d          )
+      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg_3d' , phyt_avg_3d           )
       CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max               )
+#if defined key_hadocc
       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) )
       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) )
@@ -806,8 +855,4 @@
       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         )
 #elif defined key_medusa
-      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg        )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg        )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg         )
-      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max          )
       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) )
       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd'  , trn(:,:,:,jpchd) )
@@ -826,4 +871,5 @@
       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy'  , trn(:,:,:,jpoxy) )
 #endif
+#endif
 
    END SUBROUTINE asm_bgc_bkg_wri
@@ -876,4 +922,49 @@
    !!===========================================================================
 
+   SUBROUTINE asm_bgc_unlog_3d( pbkg, pinc_log, pinc_nonlog )
+      !!------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_bgc_init_incs  ***
+      !!
+      !! ** Purpose :   convert log increments to non-log
+      !!
+      !! ** Method  :   need to account for model background,
+      !!                cannot simply do 10^log_inc. Need to:
+      !!                1) Add log_inc to log10(background) to get log10(analysis)
+      !!                2) Take 10^log10(analysis) to get analysis
+      !!                3) Subtract background from analysis to get non-log incs
+      !!
+      !! ** Action  :   return non-log increments
+      !!
+      !! References :   
+      !!------------------------------------------------------------------------
+      !!
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) :: pbkg        ! Background
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) :: pinc_log    ! Log incs
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk) :: pinc_nonlog ! Non-log incs
+      !
+      INTEGER                                         :: ji, jj, jk  ! Loop counters
+      !!
+      !!------------------------------------------------------------------------
+
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( pbkg(ji,jj,jk) > 0.0 ) THEN
+                  pinc_nonlog(ji,jj,jk) = 10**( LOG10( pbkg(ji,jj,jk) ) + &
+                     &                          pinc_log(ji,jj,jk) )      &
+                     &                    - pbkg(ji,jj,jk)
+               ELSE
+                  pinc_nonlog(ji,jj,jk) = 0.0
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+   END SUBROUTINE asm_bgc_unlog_3d
+
+   !!===========================================================================
+   !!===========================================================================
+   !!===========================================================================
+
    SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
       !!------------------------------------------------------------------------
@@ -894,23 +985,26 @@
       REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
       !
-      INTEGER                      :: jk          ! Loop counter
-      INTEGER                      :: it          ! Index
-      REAL(wp)                     :: zincwgt     ! IAU weight for current time step
-      REAL(wp)                     :: zincper     ! IAU interval in seconds
-      REAL(wp), DIMENSION(jpi,jpj) :: zmld        ! Mixed layer depth
-      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs
-      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg
-      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs
-      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg
-#if defined key_medusa
-      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
-      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
-      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs
-      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg
-      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs
-      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg
-      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs
-      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
-#endif
+      INTEGER                        :: jk          ! Loop counter
+      INTEGER                        :: it          ! Index
+      REAL(wp)                       :: zincwgt     ! IAU weight for current time step
+      REAL(wp)                       :: zincper     ! IAU interval in seconds
+      REAL(wp), DIMENSION(jpi,jpj)   :: zmld        ! Mixed layer depth
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_chltot ! Local chltot incs
+      REAL(wp), DIMENSION(jpi,jpj)   :: zbkg_chltot ! Local chltot bkg
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_phytot ! Local phytot incs
+      REAL(wp), DIMENSION(jpi,jpj)   :: zbkg_phytot ! Local phytot bkg
+#if defined key_medusa
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_chldia ! Local chldia incs
+      REAL(wp), DIMENSION(jpi,jpj)   :: zbkg_chldia ! Local chldia bkg
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_chlnon ! Local chlnon incs
+      REAL(wp), DIMENSION(jpi,jpj)   :: zbkg_chlnon ! Local chlnon bkg
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_phydia ! Local phydia incs
+      REAL(wp), DIMENSION(jpi,jpj)   :: zbkg_phydia ! Local phydia bkg
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_phynon ! Local phynon incs
+      REAL(wp), DIMENSION(jpi,jpj)   :: zbkg_phynon ! Local phynon bkg
+#endif
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zpgrow_avg_bkg ! Local pgrow_avg_bkg
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zploss_avg_bkg ! Local ploss_avg_bkg
+      REAL(wp), DIMENSION(jpi,jpj,1) :: zphyt_avg_bkg  ! Local phyt_avg_bkg
       !!------------------------------------------------------------------------
 
@@ -928,9 +1022,9 @@
             zbkg_chltot(:,:) = chl_bkg(:,:,1)
 #endif
-            CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot )
+            CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot(:,:,1) )
          ELSE IF ( ln_schltotinc ) THEN
-            zinc_chltot(:,:) = schltot_bkginc(:,:)
+            zinc_chltot(:,:,1) = schltot_bkginc(:,:)
          ELSE
-            zinc_chltot(:,:) = 0.0
+            zinc_chltot(:,:,:) = 0.0
          ENDIF
 
@@ -939,7 +1033,7 @@
          IF ( ln_slchldiainc ) THEN
             zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
-            CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia )
+            CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia(:,:,1) )
          ELSE
-            zinc_chldia(:,:) = 0.0
+            zinc_chldia(:,:,:) = 0.0
          ENDIF
 #endif
@@ -949,7 +1043,7 @@
          IF ( ln_slchlnoninc ) THEN
             zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn)
-            CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon )
+            CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon(:,:,1) )
          ELSE
-            zinc_chlnon(:,:) = 0.0
+            zinc_chlnon(:,:,:) = 0.0
          ENDIF
 #endif
@@ -962,7 +1056,7 @@
             zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
 #endif
-            CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot )
+            CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot(:,:,1) )
          ELSE
-            zinc_phytot(:,:) = 0.0
+            zinc_phytot(:,:,:) = 0.0
          ENDIF
 
@@ -971,7 +1065,7 @@
          IF ( ln_slphydiainc ) THEN
             zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd
-            CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia )
+            CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia(:,:,1) )
          ELSE
-            zinc_phydia(:,:) = 0.0
+            zinc_phydia(:,:,:) = 0.0
          ENDIF
 #endif
@@ -981,7 +1075,7 @@
          IF ( ln_slphynoninc ) THEN
             zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn
-            CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon )
+            CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon(:,:,1) )
          ELSE
-            zinc_phynon(:,:) = 0.0
+            zinc_phynon(:,:,:) = 0.0
          ENDIF
 #endif
@@ -1024,34 +1118,40 @@
 
          zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
-
-#if defined key_medusa
-         CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), &
-            &                         zinc_chltot,                         &
-            &                         ln_slchldiainc,                      &
-            &                         zinc_chldia,                         &
-            &                         ln_slchlnoninc,                      &
-            &                         zinc_chlnon,                         &
-            &                         ln_slphytotinc,                      &
-            &                         zinc_phytot,                         &
-            &                         ln_slphydiainc,                      &
-            &                         zinc_phydia,                         &
-            &                         ln_slphynoninc,                      &
-            &                         zinc_phynon,                         &
-            &                         zincper,                             &
-            &                         rn_maxchlinc, ln_phytobal, zmld,     &
-            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
-            &                         phyt_avg_bkg, mld_max_bkg,           &
-            &                         tracer_bkg, phyto2d_balinc )
+         
+         zpgrow_avg_bkg(:,:,1) = pgrow_avg_bkg(:,:)
+         zploss_avg_bkg(:,:,1) = ploss_avg_bkg(:,:)
+         zphyt_avg_bkg(:,:,1)  = phyt_avg_bkg(:,:)
+
+#if defined key_medusa
+         CALL asm_phyto_bal_medusa( 1,                                   &
+            &                       (ln_slchltotinc .OR. ln_schltotinc), &
+            &                       zinc_chltot,                         &
+            &                       ln_slchldiainc,                      &
+            &                       zinc_chldia,                         &
+            &                       ln_slchlnoninc,                      &
+            &                       zinc_chlnon,                         &
+            &                       ln_slphytotinc,                      &
+            &                       zinc_phytot,                         &
+            &                       ln_slphydiainc,                      &
+            &                       zinc_phydia,                         &
+            &                       ln_slphynoninc,                      &
+            &                       zinc_phynon,                         &
+            &                       zincper,                             &
+            &                       rn_maxchlinc, ln_phytobal, zmld,     &
+            &                       zpgrow_avg_bkg, zploss_avg_bkg,      &
+            &                       zphyt_avg_bkg, mld_max_bkg,          &
+            &                       tracer_bkg, phyto2d_balinc )
 #elif defined key_hadocc
-         CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), &
-            &                         zinc_chltot,                         &
-            &                         ln_slphytotinc,                      &
-            &                         zinc_phytot,                         &
-            &                         zincper,                             &
-            &                         rn_maxchlinc, ln_phytobal, zmld,     &
-            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
-            &                         phyt_avg_bkg, mld_max_bkg,           &
-            &                         cchl_p_bkg(:,:,1),                   &
-            &                         tracer_bkg, phyto2d_balinc )
+         CALL asm_phyto_bal_hadocc( 1,                                   &
+            &                       (ln_slchltotinc .OR. ln_schltotinc), &
+            &                       zinc_chltot,                         &
+            &                       ln_slphytotinc,                      &
+            &                       zinc_phytot,                         &
+            &                       zincper,                             &
+            &                       rn_maxchlinc, ln_phytobal, zmld,     &
+            &                       zpgrow_avg_bkg, zploss_avg_bkg,      &
+            &                       zphyt_avg_bkg, mld_max_bkg,          &
+            &                       cchl_p_bkg(:,:,1),                   &
+            &                       tracer_bkg, phyto2d_balinc )
 #else
          CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
@@ -1166,12 +1266,10 @@
       INTEGER                          :: ji, jj, jk   ! Loop counters
       INTEGER                          :: it           ! Index
+      REAL(wp)                         :: zincper      ! IAU interval in seconds
       REAL(wp)                         :: zincwgt      ! IAU weight for timestep
-      REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn
-      REAL(wp)                         :: zfrac_chd    ! Fraction of jpchd
-      REAL(wp)                         :: zrat_phn_chn ! jpphn:jpchn ratio
-      REAL(wp)                         :: zrat_phd_chd ! jpphd:jpchd ratio
-      REAL(wp)                         :: zrat_pds_chd ! jppds:jpchd ratio
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc      ! Chlorophyll increments
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl      ! Chlorophyll background
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zinc_chltot  ! Chlorophyll increments
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbkg_chltot  ! Chlorophyll background
+      REAL(wp), DIMENSION(jpi,jpj)     :: zdummy_2d    ! Dummy array for call
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdummy_3d    ! Dummy array for call
       !!------------------------------------------------------------------------
 
@@ -1179,68 +1277,49 @@
 
          IF ( ln_plchltotinc ) THEN
-            ! Convert log10(chlorophyll) increment back to a chlorophyll increment
-            ! In order to transform logchl incs to chl incs, need to account for model
-            ! background, cannot simply do 10^logchl_bkginc. Need to:
-            ! 1) Add logchl inc to log10(background) to get log10(analysis)
-            ! 2) Take 10^log10(analysis) to get analysis
-            ! 3) Subtract background from analysis to get chl incs
-            ! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
-#if defined key_medusa
-            bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
+#if defined key_medusa
+            zbkg_chltot(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
 #elif defined key_hadocc
-            bkg_chl(:,:,:) = chl_bkg(:,:,:)
-#endif
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
-                        chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
-                        IF ( rn_maxchlinc > 0.0 ) THEN
-                           chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
-                        ENDIF
-                     ELSE
-                        chl_inc(ji,jj,jk) = 0.0
-                     ENDIF
-                  END DO
-               END DO
-            END DO
+            zbkg_chltot(:,:,:) = chl_bkg(:,:,:)
+#endif
+            CALL asm_bgc_unlog_3d( zbkg_chltot, plchltot_bkginc, zinc_chltot )
          ELSE IF ( ln_pchltotinc ) THEN
-            DO jk = 1, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     IF ( rn_maxchlinc > 0.0 ) THEN
-                        chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
-                     ELSE
-                        chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
-                     ENDIF
-                  END DO
-               END DO
-            END DO
-         ENDIF
-
-#if defined key_medusa
-         ! Loop over each grid point partioning the increments based on existing ratios
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
-                     zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
-                     zfrac_chd = 1.0 - zfrac_chn
-                     phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
-                     phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
-                     zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
-                     zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
-                     zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
-                     phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
-                     phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
-                     phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
-                  ENDIF
-               END DO
-            END DO
-         END DO
+            zinc_chltot(:,:,:) = pchltot_bkginc(:,:,:)
+         ENDIF
+
+         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
+
+#if defined key_medusa
+         CALL asm_phyto_bal_medusa( jpk,                                  &
+            &                       (ln_plchltotinc .OR. ln_pchltotinc),  &
+            &                       zinc_chltot,                          &
+            &                       .FALSE.,                              &
+            &                       zdummy_3d,                            &
+            &                       .FALSE.,                              &
+            &                       zdummy_3d,                            &
+            &                       .FALSE.,                              &
+            &                       zdummy_3d,                            &
+            &                       .FALSE.,                              &
+            &                       zdummy_3d,                            &
+            &                       .FALSE.,                              &
+            &                       zdummy_3d,                            &
+            &                       zincper,                              &
+            &                       rn_maxchlinc, ln_phytobal, zdummy_2d, &
+            &                       pgrow_avg_3d_bkg, ploss_avg_3d_bkg,   &
+            &                       phyt_avg_3d_bkg, mld_max_bkg,         &
+            &                       tracer_bkg, phyto3d_balinc )
 #elif defined key_hadocc
-         phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
-#else
-         CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
+         CALL asm_phyto_bal_hadocc( jpk,                                  &
+            &                       (ln_plchltotinc .OR. ln_pchltotinc),  &
+            &                       zinc_chltot,                          &
+            &                       .FALSE.,                              &
+            &                       zdummy_3d,                            &
+            &                       zincper,                              &
+            &                       rn_maxchlinc, ln_phytobal, zdummy_2d, &
+            &                       pgrow_avg_3d_bkg, ploss_avg_3d_bkg,   &
+            &                       phyt_avg_3d_bkg, mld_max_bkg,         &
+            &                       cchl_p_bkg,                           &
+            &                       tracer_bkg, phyto3d_balinc )
+#else
+         CALL ctl_stop( 'Attempting to assimilate phyto3d data, ', &
             &           'but not defined a biogeochemical model' )
 #endif
@@ -1423,6 +1502,12 @@
          ! Account for phytoplankton balancing if required
          IF ( ln_phytobal ) THEN
-            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
-            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
+            IF ( ALLOCATED(phyto2d_balinc) ) THEN
+               dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
+               alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
+            ENDIF
+            IF ( ALLOCATED(phyto3d_balinc) ) THEN
+               dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto3d_balinc(:,:,1,jpdic)
+               alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto3d_balinc(:,:,1,jpalk)
+            ENDIF
          ELSE
             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
@@ -1437,6 +1522,12 @@
          ! Account for phytoplankton balancing if required
          IF ( ln_phytobal ) THEN
-            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
-            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
+            IF ( ALLOCATED(phyto2d_balinc) ) THEN
+               dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
+               alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
+            ENDIF
+            IF ( ALLOCATED(phyto3d_balinc) ) THEN
+               dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto3d_balinc(:,:,1,jp_had_dic)
+               alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto3d_balinc(:,:,1,jp_had_alk)
+            ENDIF
          ELSE
             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
@@ -1677,8 +1768,16 @@
          ! Account for phytoplankton balancing if required
          IF ( ln_phytobal ) THEN
-            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
-            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
-            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
-            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
+            IF ( ALLOCATED(phyto2d_balinc) ) THEN
+               dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
+               alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
+               din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
+               sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
+            ENDIF
+            IF ( ALLOCATED(phyto3d_balinc) ) THEN
+               dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto3d_balinc(:,:,:,jpdic)
+               alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto3d_balinc(:,:,:,jpalk)
+               din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto3d_balinc(:,:,:,jpdin)
+               sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto3d_balinc(:,:,:,jpsil)
+            ENDIF
          ELSE
             dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
@@ -1871,8 +1970,8 @@
             it = jpdin
 #endif
-            IF ( ln_phytobal ) THEN
+            IF ( ALLOCATED(phyto2d_balinc) ) THEN
                pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
-            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+            IF ( ALLOCATED(phyto3d_balinc) ) THEN
                pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
             ENDIF
@@ -1891,8 +1990,8 @@
 #if defined key_medusa
             it = jpsil
-            IF ( ln_phytobal ) THEN
+            IF ( ALLOCATED(phyto2d_balinc) ) THEN
                psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
-            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+            IF ( ALLOCATED(phyto3d_balinc) ) THEN
                psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
             ENDIF
@@ -1915,8 +2014,8 @@
             it = jpdic
 #endif
-            IF ( ln_phytobal ) THEN
+            IF ( ALLOCATED(phyto2d_balinc) ) THEN
                pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
-            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+            IF ( ALLOCATED(phyto3d_balinc) ) THEN
                pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
             ENDIF
@@ -1939,8 +2038,8 @@
             it = jpalk
 #endif
-            IF ( ln_phytobal ) THEN
+            IF ( ALLOCATED(phyto2d_balinc) ) THEN
                palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
-            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+            IF ( ALLOCATED(phyto3d_balinc) ) THEN
                palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
             ENDIF
@@ -1959,8 +2058,8 @@
 #if defined key_medusa
             it = jpoxy
-            IF ( ln_phytobal ) THEN
+            IF ( ALLOCATED(phyto2d_balinc) ) THEN
                po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
-            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
+            IF ( ALLOCATED(phyto3d_balinc) ) THEN
                po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
             ENDIF
