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 13353)
+++ /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
Index: anches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_hadocc.F90
===================================================================
--- /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_hadocc.F90	(revision 13353)
+++ 	(revision )
@@ -1,291 +1,0 @@
-MODULE asmphyto2dbal_hadocc
-   !!======================================================================
-   !!                       ***  MODULE asmphyto2dbal_hadocc  ***
-   !! Calculate increments to HadOCC based on surface phyto2d increments
-   !!
-   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
-   !! For licensing reasons this is kept in its own internal Met Office
-   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
-   !! repository, and must be merged in when building.
-   !!
-   !!======================================================================
-   !! History :  3.6  ! 2017-08  (D. Ford)     Adapted from bioanal.F90
-   !!----------------------------------------------------------------------
-#if defined key_asminc && defined key_hadocc
-   !!----------------------------------------------------------------------
-   !! 'key_asminc'          : assimilation increment interface
-   !! 'key_hadocc'          : HadOCC model
-   !!----------------------------------------------------------------------
-   !! asm_phyto2d_bal_hadocc : routine to calculate increments to HadOCC
-   !!----------------------------------------------------------------------
-   USE par_kind,      ONLY: wp             ! kind parameters
-   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
-   USE dom_oce,       ONLY: gdepw_n        ! domain information
-   USE iom                                 ! i/o
-   USE par_hadocc                          ! HadOCC parameters
-   USE had_bgc_stnd,  ONLY: kmt            ! HadOCC parameters
-   USE had_bgc_const                       ! HadOCC parameters
-   USE par_trc,       ONLY: jptra          ! Tracer parameters
-   USE bioanalysis                         ! Nitrogen balancing
-
-   IMPLICIT NONE
-   PRIVATE                   
-
-   PUBLIC asm_phyto2d_bal_hadocc
-
-   ! Default values for biological assimilation parameters
-   ! Should match Hemmings et al. (2008)
-   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
-   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
-   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
-   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
-   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
-   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
-   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
-   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
-   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
-   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
-   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
-   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
-   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
-   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
-   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
-   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
-   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
-   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
-   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
-   !
-   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
-   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
-   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
-   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
-   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
-   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
-   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
-   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
-
-CONTAINS
-
-   SUBROUTINE asm_phyto2d_bal_hadocc( ld_chltot,                      &
-      &                               pinc_chltot,                    &
-      &                               ld_phytot,                      &
-      &                               pinc_phytot,                    &
-      &                               pincper,                        &
-      &                               p_maxchlinc, ld_phytobal, pmld, &
-      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
-      &                               phyt_avg_bkg, mld_max_bkg,      &
-      &                               cchl_p_bkg,                     &
-      &                               tracer_bkg, phyto2d_balinc )
-      !!---------------------------------------------------------------------------
-      !!                    ***  ROUTINE asm_phyto2d_bal_hadocc  ***
-      !!
-      !! ** Purpose :   calculate increments to HadOCC from 2d phytoplankton increments
-      !!
-      !! ** Method  :   call nitrogen balancing scheme
-      !!
-      !! ** Action  :   populate phyto2d_balinc
-      !!
-      !! References :   Hemmings et al., 2008, J. Mar. Res.
-      !!                Ford et al., 2012, Ocean Sci.
-      !!---------------------------------------------------------------------------
-      !!
-      LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments
-      LOGICAL,  INTENT(in   )                               :: ld_phytot      ! Assim phytot y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phytot    ! phytot increments
-      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
-      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
-      LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: cchl_p_bkg     ! Surface C:Chl
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
-      !!
-      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
-      INTEGER                                               :: jkmax          ! Loop index
-      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
-      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
-      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
-      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
-      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
-      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
-      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
-      !!---------------------------------------------------------------------------
-
-      IF ( ( .NOT. ld_chltot ) .AND. ( .NOT. ld_phytot ) ) THEN
-         CALL ctl_stop( ' Trying to do phyto2d balancing but nothing to assimilate' )
-      ENDIF
-      
-      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
-      IF ( p_maxchlinc > 0.0 ) THEN
-         IF ( ld_chltot ) THEN
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
-               END DO
-            END DO
-         ENDIF
-      ENDIF
-
-      IF ( ld_phytot ) THEN
-         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
-      ENDIF
-      
-      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
-
-         ! Set up model parameters to be passed into Hemmings balancing routine
-         modparm(1)  = grow_sat
-         modparm(2)  = psmax
-         modparm(3)  = par
-         modparm(4)  = alpha
-         modparm(5)  = resp_rate
-         modparm(6)  = pmort_rate
-         modparm(7)  = phyto_min
-         modparm(8)  = z_mort_1
-         modparm(9)  = z_mort_2
-         modparm(10) = c2n_p
-         modparm(11) = c2n_z
-         modparm(12) = c2n_d
-         modparm(13) = graze_threshold
-         modparm(14) = holling_coef
-         modparm(15) = graze_sat
-         modparm(16) = graze_max
-
-         ! Set up assimilation parameters to be passed into balancing routine
-         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
-         assimparm(2)  = balnutext
-         assimparm(3)  = balnutmin
-         assimparm(4)  = r
-         assimparm(5)  = beta_g
-         assimparm(6)  = beta_l
-         assimparm(7)  = beta_m
-         assimparm(8)  = a_g
-         assimparm(9)  = a_l
-         assimparm(10) = a_m
-         assimparm(11) = zfracb0
-         assimparm(12) = zfracb1
-         assimparm(13) = qrfmax
-         assimparm(14) = qafmax
-         assimparm(15) = zrfmax
-         assimparm(16) = zafmax
-         assimparm(17) = prfmax
-         assimparm(18) = incphymin
-         assimparm(19) = integnstep
-         assimparm(20) = pthreshold
-
-         ! Set up external tracer indices array bstate
-         i_tracer(1) = 1   ! nutrient
-         i_tracer(2) = 2   ! phytoplankton
-         i_tracer(3) = 3   ! zooplankton
-         i_tracer(4) = 4   ! detritus
-         i_tracer(5) = 5   ! DIC
-         i_tracer(6) = 6   ! Alkalinity
-
-         ! Set background state
-         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_had_nut)
-         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_had_phy)
-         bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jp_had_zoo)
-         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_had_pdn)
-         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_had_dic)
-         bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jp_had_alk)
-
-         ! Call nitrogen balancing routine
-         CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,               &
-            &               n2be_p, n2be_z, n2be_d, assimparm,                          &
-            &               INT(pincper), 1, kmt(:,:), tmask(:,:,:),                    &
-            &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p_bkg(:,:), &
-            &               nbal_active, phyt_avg_bkg(:,:),                             &
-            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),          &
-            &               subsurf_active, deepneg_active,                             &
-            &               deeppos_active, nutprof_active,                             &
-            &               bstate, outincs,                                            &
-            &               diag_active, diag,                                          &
-            &               diag_fulldepth_active, diag_fulldepth )
-
-         ! Save balancing increments
-         phyto2d_balinc(:,:,:,jp_had_nut) = outincs(:,:,:,i_tracer(1))
-         phyto2d_balinc(:,:,:,jp_had_phy) = outincs(:,:,:,i_tracer(2))
-         phyto2d_balinc(:,:,:,jp_had_zoo) = outincs(:,:,:,i_tracer(3))
-         phyto2d_balinc(:,:,:,jp_had_pdn) = outincs(:,:,:,i_tracer(4))
-         phyto2d_balinc(:,:,:,jp_had_dic) = outincs(:,:,:,i_tracer(5))
-         phyto2d_balinc(:,:,:,jp_had_alk) = outincs(:,:,:,i_tracer(6))
-      
-      ELSE   ! No nitrogen balancing
-      
-         ! Initialise phytoplankton increment to zero
-         phyto2d_balinc(:,:,:,jp_had_phy) = 0.0
-         
-         ! Convert surface chlorophyll increment to phytoplankton nitrogen
-         phyto2d_balinc(:,:,1,jp_had_phy) = ( cchl_p_bkg(:,:) / (mw_carbon * c2n_p) ) * pinc_chltot(:,:)
-         
-         ! Propagate through mixed layer
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               jkmax = jpk-1
-               DO jk = jpk-1, 1, -1
-                  IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
-                     & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
-                     pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
-                     jkmax = jk
-                  ENDIF
-               END DO
-               !
-               DO jk = 2, jkmax
-                  phyto2d_balinc(ji,jj,jk,jp_had_phy) = phyto2d_balinc(ji,jj,1,jp_had_phy)
-               END DO
-               !
-            END DO
-         END DO
-
-         ! Set other balancing increments to zero
-         phyto2d_balinc(:,:,:,jp_had_nut) = 0.0
-         phyto2d_balinc(:,:,:,jp_had_zoo) = 0.0
-         phyto2d_balinc(:,:,:,jp_had_pdn) = 0.0
-         phyto2d_balinc(:,:,:,jp_had_dic) = 0.0
-         phyto2d_balinc(:,:,:,jp_had_alk) = 0.0
-      
-      ENDIF
-
-   END SUBROUTINE asm_phyto2d_bal_hadocc
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option : Empty routine
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE asm_phyto2d_bal_hadocc( ld_chltot,                      &
-      &                               pinc_chltot,                    &
-      &                               ld_phytot,                      &
-      &                               pinc_phytot,                    &
-      &                               pincper,                        &
-      &                               p_maxchlinc, ld_phytobal, pmld, &
-      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
-      &                               phyt_avg_bkg, mld_max_bkg,      &
-      &                               cchl_p_bkg,                     &
-      &                               tracer_bkg, phyto2d_balinc )
-      LOGICAL :: ld_chltot
-      REAL    :: pinc_chltot(:,:)
-      LOGICAL :: ld_phytot
-      REAL    :: pinc_phytot(:,:)
-      REAL    :: pincper
-      REAL    :: p_maxchlinc
-      LOGICAL :: ld_phytobal
-      REAL    :: pmld(:,:)
-      REAL    :: pgrow_avg_bkg(:,:)
-      REAL    :: ploss_avg_bkg(:,:)
-      REAL    :: phyt_avg_bkg(:,:)
-      REAL    :: mld_max_bkg(:,:)
-      REAL    :: cchl_p_bkg(:,:)
-      REAL    :: tracer_bkg(:,:,:,:)
-      REAL    :: phyto2d_balinc(:,:,:,:)
-      WRITE(*,*) 'asm_phyto2d_bal_hadocc: You should not have seen this print! error?'
-   END SUBROUTINE asm_phyto2d_bal_hadocc
-#endif
-
-   !!======================================================================
-END MODULE asmphyto2dbal_hadocc
Index: anches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90
===================================================================
--- /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90	(revision 13353)
+++ 	(revision )
@@ -1,510 +1,0 @@
-MODULE asmphyto2dbal_medusa
-   !!======================================================================
-   !!                       ***  MODULE asmphyto2dbal_medusa  ***
-   !! Calculate increments to MEDUSA based on surface phyto2d increments
-   !!
-   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
-   !! For licensing reasons this is kept in its own internal Met Office
-   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
-   !! repository, and must be merged in when building.
-   !!
-   !!======================================================================
-   !! History :  3.6  ! 2017-08 (D. Ford)  Adapted from asmphyto2dbal_hadocc
-   !!----------------------------------------------------------------------
-#if defined key_asminc && defined key_medusa
-   !!----------------------------------------------------------------------
-   !! 'key_asminc'          : assimilation increment interface
-   !! 'key_medusa'          : MEDUSA model
-   !!----------------------------------------------------------------------
-   !! asm_phyto2d_bal_medusa : routine to calculate increments to MEDUSA
-   !!----------------------------------------------------------------------
-   USE par_kind,      ONLY: wp             ! kind parameters
-   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
-   USE dom_oce,       ONLY: gdepw_n        ! domain information
-   USE zdftmx,        ONLY: ln_tmx_itf, &  ! Indonesian Throughflow
-      &                     mask_itf       ! tidal mixing mask
-   USE iom                                 ! i/o
-   USE sms_medusa                          ! MEDUSA parameters
-   USE par_medusa                          ! MEDUSA parameters
-   USE par_trc,       ONLY: jptra          ! Tracer parameters
-   USE bioanalysis                         ! Nitrogen balancing
-
-   IMPLICIT NONE
-   PRIVATE                   
-
-   PUBLIC asm_phyto2d_bal_medusa
-
-   ! Default values for biological assimilation parameters
-   ! Should match Hemmings et al. (2008)
-   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
-   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
-   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
-   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
-   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
-   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
-   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
-   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
-   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
-   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
-   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
-   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
-   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
-   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
-   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
-   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
-   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
-   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
-   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
-   !
-   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
-   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
-   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
-   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
-   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
-   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
-   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
-   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
-
-CONTAINS
-
-   SUBROUTINE asm_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 )
-      !!---------------------------------------------------------------------------
-      !!                    ***  ROUTINE asm_phyto2d_bal_medusa  ***
-      !!
-      !! ** Purpose :   calculate increments to MEDUSA from 2d phytoplankton increments
-      !!
-      !! ** Method  :   average up MEDUSA to look like HadOCC
-      !!                call nitrogen balancing scheme
-      !!                separate back out to MEDUSA
-      !!
-      !! ** Action  :   populate phyto2d_balinc
-      !!
-      !! References :   Hemmings et al., 2008, J. Mar. Res.
-      !!                Ford et al., 2012, Ocean Sci.
-      !!---------------------------------------------------------------------------
-      !!
-      LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments
-      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
-      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
-      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
-      LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
-      !!
-      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
-      INTEGER                                               :: jkmax          ! Loop index
-      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
-      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
-      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
-      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
-      REAL(wp)                                              :: zfrac          ! 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),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
-      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
-      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
-      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
-      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
-      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
-      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
-      !!---------------------------------------------------------------------------
-
-      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
-      IF ( p_maxchlinc > 0.0 ) THEN
-         IF ( ld_chltot ) THEN
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
-               END DO
-            END DO
-         ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
-            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)
-               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
-      
-      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
-
-         ! Set up model parameters to be passed into Hemmings balancing routine.
-         ! For now these are hardwired to the standard HadOCC parameter values
-         ! (except C:N ratios) as this is what the scheme was developed for.
-         ! Obviously, HadOCC and MEDUSA are rather different models, so this
-         ! isn't ideal, but there's not always direct analogues between the two
-         ! parameter sets, so it's the easiest way to get something running.
-         ! In the longer term, some serious MarMOT-based development is required.
-         modparm(1)  = 0.1                               ! grow_sat
-         modparm(2)  = 2.0                               ! psmax
-         modparm(3)  = 0.845                             ! par
-         modparm(4)  = 0.02                              ! alpha
-         modparm(5)  = 0.05                              ! resp_rate
-         modparm(6)  = 0.05                              ! pmort_rate
-         modparm(7)  = 0.01                              ! phyto_min
-         modparm(8)  = 0.05                              ! z_mort_1
-         modparm(9)  = 1.0                               ! z_mort_2
-         modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p
-         modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z
-         modparm(12) = xthetad                           ! c2n_d
-         modparm(13) = 0.01                              ! graze_threshold
-         modparm(14) = 2.0                               ! holling_coef
-         modparm(15) = 0.5                               ! graze_sat
-         modparm(16) = 2.0                               ! graze_max
-
-         ! Set up assimilation parameters to be passed into balancing routine
-         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
-         assimparm(2)  = balnutext
-         assimparm(3)  = balnutmin
-         assimparm(4)  = r
-         assimparm(5)  = beta_g
-         assimparm(6)  = beta_l
-         assimparm(7)  = beta_m
-         assimparm(8)  = a_g
-         assimparm(9)  = a_l
-         assimparm(10) = a_m
-         assimparm(11) = zfracb0
-         assimparm(12) = zfracb1
-         assimparm(13) = qrfmax
-         assimparm(14) = qafmax
-         assimparm(15) = zrfmax
-         assimparm(16) = zafmax
-         assimparm(17) = prfmax
-         assimparm(18) = incphymin
-         assimparm(19) = integnstep
-         assimparm(20) = pthreshold
-
-         ! Set up external tracer indices array bstate
-         i_tracer(1) = 1   ! nutrient
-         i_tracer(2) = 2   ! phytoplankton
-         i_tracer(3) = 3   ! zooplankton
-         i_tracer(4) = 4   ! detritus
-         i_tracer(5) = 5   ! DIC
-         i_tracer(6) = 6   ! Alkalinity
-
-         ! Set background state
-         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin)
-         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd)
-         bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme)
-         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet)
-         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic)
-         bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk)
-
-         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
-         ! and nitrogen to biomass equivalent for PZD
-         ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA
-         cchl_p(:,:) = 0.0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN
-                  cchl_p(ji,jj) = xmassc * ( ( tracer_bkg(ji,jj,1,jpphn) * xthetapn ) +      &
-                     &                       ( tracer_bkg(ji,jj,1,jpphd) * xthetapd )   ) /  &
-                     &            ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) )
-               ENDIF
-            END DO
-         END DO
-         n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) )
-         n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) )
-         n2be_d = 14.01 + ( xmassc * xthetad )
-
-         ! Call nitrogen balancing routine
-         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
-            &               n2be_p, n2be_z, n2be_d, assimparm,                      &
-            &               INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
-            &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p(:,:), &
-            &               nbal_active, phyt_avg_bkg(:,:),                         &
-            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      &
-            &               subsurf_active, deepneg_active,                         &
-            &               deeppos_active, nutprof_active,                         &
-            &               bstate, outincs,                                        &
-            &               diag_active, diag,                                      &
-            &               diag_fulldepth_active, diag_fulldepth )
-         
-         ! Loop over each grid point partioning the increments
-         phyto2d_balinc(:,:,:,:) = 0.0
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-
-                  ! Phytoplankton
-                  IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. &
-                     & ( tracer_bkg(ji,jj,jk,jpphd) > 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
-                        ! 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)
-                     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)
-                     
-                     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
-                  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))
-
-                  ! DIC straight from balancing scheme
-                  phyto2d_balinc(ji,jj,jk,jpdic) = 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
-                  
-               END DO
-            END DO
-         END DO
-      
-      ELSE   ! No nitrogen balancing
-      
-         ! Initialise individual chlorophyll increments to zero
-         phyto2d_balinc(:,:,:,jpchn) = 0.0
-         phyto2d_balinc(:,:,:,jpchd) = 0.0
-         
-         ! 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 ( 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
-                  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)
-                  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
-                  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
-                  ENDIF
-               ENDIF
-            END DO
-         END DO
-         
-         ! Propagate through mixed layer
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               jkmax = jpk-1
-               DO jk = jpk-1, 1, -1
-                  IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
-                     & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
-                     pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
-                     jkmax = jk
-                  ENDIF
-               END DO
-               !
-               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)
-               END DO
-               !
-            END DO
-         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
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option : Empty routine
-   !!----------------------------------------------------------------------
-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 )
-      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(:,:)
-      REAL    :: pincper
-      REAL    :: p_maxchlinc
-      LOGICAL :: ld_phytobal
-      REAL    :: pmld(:,:)
-      REAL    :: pgrow_avg_bkg(:,:)
-      REAL    :: ploss_avg_bkg(:,:)
-      REAL    :: phyt_avg_bkg(:,:)
-      REAL    :: mld_max_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
-#endif
-
-   !!======================================================================
-END MODULE asmphyto2dbal_medusa
Index: /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_hadocc.F90
===================================================================
--- /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_hadocc.F90	(revision 13354)
+++ /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_hadocc.F90	(revision 13354)
@@ -0,0 +1,346 @@
+MODULE asmphytobal_hadocc
+   !!======================================================================
+   !!                       ***  MODULE asmphytobal_hadocc  ***
+   !! Calculate increments to HadOCC based on surface phyto increments
+   !!
+   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
+   !! For licensing reasons this is kept in its own internal Met Office
+   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
+   !! repository, and must be merged in when building.
+   !!
+   !!======================================================================
+   !! History :  3.6  ! 2017-08  (D. Ford)     Adapted from bioanal.F90
+   !!----------------------------------------------------------------------
+#if defined key_asminc && defined key_hadocc
+   !!----------------------------------------------------------------------
+   !! 'key_asminc'          : assimilation increment interface
+   !! 'key_hadocc'          : HadOCC model
+   !!----------------------------------------------------------------------
+   !! asm_phyto_bal_hadocc : routine to calculate increments to HadOCC
+   !!----------------------------------------------------------------------
+   USE par_kind,      ONLY: wp             ! kind parameters
+   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
+   USE dom_oce,       ONLY: gdepw_n        ! domain information
+   USE iom                                 ! i/o
+   USE par_hadocc                          ! HadOCC parameters
+   USE had_bgc_stnd,  ONLY: kmt            ! HadOCC parameters
+   USE had_bgc_const                       ! HadOCC parameters
+   USE par_trc,       ONLY: jptra          ! Tracer parameters
+   USE bioanalysis                         ! Nitrogen balancing
+
+   IMPLICIT NONE
+   PRIVATE                   
+
+   PUBLIC asm_phyto_bal_hadocc
+
+   ! Default values for biological assimilation parameters
+   ! Should match Hemmings et al. (2008)
+   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
+   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
+   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
+   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
+   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
+   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
+   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
+   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
+   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
+   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
+   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
+   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
+   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
+   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
+   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
+   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
+   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
+   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
+   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
+   !
+   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
+   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
+   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
+   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
+   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
+   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
+   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
+   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
+
+CONTAINS
+
+   SUBROUTINE asm_phyto_bal_hadocc( kdeps,                              &
+      &                             ld_chltot,                          &
+      &                             pinc_chltot_3d,                     &
+      &                             ld_phytot,                          &
+      &                             pinc_phytot_3d,                     &
+      &                             pincper,                            &
+      &                             p_maxchlinc, ld_phytobal, pmld,     &
+      &                             pgrow_avg_bkg_3d, ploss_avg_bkg_3d, &
+      &                             phyt_avg_bkg_3d, mld_max_bkg,       &
+      &                             cchl_p_bkg_3d,                      &
+      &                             tracer_bkg, phyto_balinc )
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_phyto_bal_hadocc  ***
+      !!
+      !! ** Purpose :   calculate increments to HadOCC from 2d phytoplankton increments
+      !!
+      !! ** Method  :   call nitrogen balancing scheme
+      !!
+      !! ** Action  :   populate phyto_balinc
+      !!
+      !! References :   Hemmings et al., 2008, J. Mar. Res.
+      !!                Ford et al., 2012, Ocean Sci.
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER,  INTENT(in   )                               :: kdeps            ! No. inc deps 1 or jpk
+      LOGICAL,  INTENT(in   )                               :: ld_chltot        ! Assim chltot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_chltot_3d   ! chltot increments
+      LOGICAL,  INTENT(in   )                               :: ld_phytot        ! Assim phytot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_phytot_3d   ! phytot increments
+      REAL(wp), INTENT(in   )                               :: pincper          ! Assimilation period
+      REAL(wp), INTENT(in   )                               :: p_maxchlinc      ! Max chl increment
+      LOGICAL,  INTENT(in   )                               :: ld_phytobal      ! Balancing y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld             ! Mixed layer depth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: pgrow_avg_bkg_3d ! Avg phyto growth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: ploss_avg_bkg_3d ! Avg phyto loss
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: phyt_avg_bkg_3d  ! Avg phyto
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg      ! Max MLD
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: cchl_p_bkg_3d    ! C:Chl
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg       ! State variables
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto_balinc     ! Balancing increments
+      !!
+      INTEGER                                               :: ji, jj, jk, jn    ! Loop counters
+      INTEGER                                               :: jkmax             ! Loop index
+      INTEGER,                 DIMENSION(6)                 :: i_tracer          ! Tracer indices
+      REAL(wp),                DIMENSION(16)                :: modparm           ! Model parameters
+      REAL(wp),                DIMENSION(20)                :: assimparm         ! Assimilation parameters
+      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: bstate_2d         ! Background state (2D)
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate_3d         ! Background state (3D)
+      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: outincs_2d        ! Balancing increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs_3d        ! Balancing increments (3D)
+      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag              ! Depth-indep diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj,1,22)      :: diag_fulldepth_2d ! Full-depth diagnostics (2D)
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth_3d ! Full-depth diagnostics (3D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p_bkg_2d     ! C:Chl for total phy (2D)
+      REAL(wp),                DIMENSION(jpi,jpj,1)         :: tmask_2d          ! Single-level tmask
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_2d    ! chltot increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_phytot_2d    ! phytot increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pgrow_avg_bkg_2d  ! Avg phyto growth (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: ploss_avg_bkg_2d  ! Avg phyto loss (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: phyt_avg_bkg_2d   ! Avg phyto (2D)
+      !!---------------------------------------------------------------------------
+
+      IF ( ( .NOT. ld_chltot ) .AND. ( .NOT. ld_phytot ) ) THEN
+         CALL ctl_stop( ' Trying to do phyto balancing but nothing to assimilate' )
+      ENDIF
+      
+      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      IF ( p_maxchlinc > 0.0 ) THEN
+         IF ( ld_chltot ) THEN
+            DO jk = 1, kdeps
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     pinc_chltot_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_3d(ji,jj,jk), p_maxchlinc ) )
+                  END DO
+               END DO
+            END DO
+         ENDIF
+      ENDIF
+
+      IF ( ld_phytot ) THEN
+         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
+      ENDIF
+      
+      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
+
+         ! Set up model parameters to be passed into Hemmings balancing routine
+         modparm(1)  = grow_sat
+         modparm(2)  = psmax
+         modparm(3)  = par
+         modparm(4)  = alpha
+         modparm(5)  = resp_rate
+         modparm(6)  = pmort_rate
+         modparm(7)  = phyto_min
+         modparm(8)  = z_mort_1
+         modparm(9)  = z_mort_2
+         modparm(10) = c2n_p
+         modparm(11) = c2n_z
+         modparm(12) = c2n_d
+         modparm(13) = graze_threshold
+         modparm(14) = holling_coef
+         modparm(15) = graze_sat
+         modparm(16) = graze_max
+
+         ! Set up assimilation parameters to be passed into balancing routine
+         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
+         assimparm(2)  = balnutext
+         assimparm(3)  = balnutmin
+         assimparm(4)  = r
+         assimparm(5)  = beta_g
+         assimparm(6)  = beta_l
+         assimparm(7)  = beta_m
+         assimparm(8)  = a_g
+         assimparm(9)  = a_l
+         assimparm(10) = a_m
+         assimparm(11) = zfracb0
+         assimparm(12) = zfracb1
+         assimparm(13) = qrfmax
+         assimparm(14) = qafmax
+         assimparm(15) = zrfmax
+         assimparm(16) = zafmax
+         assimparm(17) = prfmax
+         assimparm(18) = incphymin
+         assimparm(19) = integnstep
+         assimparm(20) = pthreshold
+
+         ! Set up external tracer indices array bstate
+         i_tracer(1) = 1   ! nutrient
+         i_tracer(2) = 2   ! phytoplankton
+         i_tracer(3) = 3   ! zooplankton
+         i_tracer(4) = 4   ! detritus
+         i_tracer(5) = 5   ! DIC
+         i_tracer(6) = 6   ! Alkalinity
+
+         ! Set background state
+         bstate_3d(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_had_nut)
+         bstate_3d(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_had_phy)
+         bstate_3d(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jp_had_zoo)
+         bstate_3d(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_had_pdn)
+         bstate_3d(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_had_dic)
+         bstate_3d(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jp_had_alk)
+
+         ! Call nitrogen balancing routine
+         IF (kdeps == 1) THEN
+            pinc_chltot_2d(:,:)   = pinc_chltot_3d(:,:,1)
+            cchl_p_bkg_2d(:,:)    = cchl_p_bkg_3d(:,:,1)
+            phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg_3d(:,:,1)
+            pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg_3d(:,:,1)
+            ploss_avg_bkg_2d(:,:) = ploss_avg_bkg_3d(:,:,1)
+            
+            CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,               &
+               &               n2be_p, n2be_z, n2be_d, assimparm,                          &
+               &               INT(pincper), 1, kmt(:,:), tmask(:,:,:),                    &
+               &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot_2d(:,:), cchl_p_bkg_2d(:,:), &
+               &               nbal_active, phyt_avg_bkg_2d(:,:),                          &
+               &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:),    &
+               &               subsurf_active, deepneg_active,                             &
+               &               deeppos_active, nutprof_active,                             &
+               &               bstate_3d, outincs_3d,                                      &
+               &               diag_active, diag,                                          &
+               &               diag_fulldepth_active, diag_fulldepth_3d )
+         ELSE
+            pmld(:,:) = 0.5
+            
+            DO jk = 1, kdeps
+               pinc_chltot_2d(:,:)   = pinc_chltot_3d(:,:,jk)
+               cchl_p_bkg_2d(:,:)    = cchl_p_bkg_3d(:,:,jk)
+               phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg_3d(:,:,jk)
+               pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg_3d(:,:,jk)
+               ploss_avg_bkg_2d(:,:) = ploss_avg_bkg_3d(:,:,jk)
+               tmask_2d(:,:,1)       = tmask(:,:,jk)
+               bstate_2d(:,:,1,:)    = bstate_3d(:,:,jk,:)
+               outincs_2d(:,:,:,:)   = 0.0
+
+               CALL bio_analysis( jpi, jpj, 1, gdepw_n(:,:,2), i_tracer, modparm,            &
+                  &               n2be_p, n2be_z, n2be_d, assimparm,                         &
+                  &               INT(pincper), 1, INT(SUM(tmask_2d,3)), tmask_2d(:,:,:),    &
+                  &               pmld(:,:), pmld(:,:), pinc_chltot_2d(:,:), cchl_p_bkg_2d(:,:), &
+                  &               nbal_active, phyt_avg_bkg_2d(:,:),                         &
+                  &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:),   &
+                  &               subsurf_active, deepneg_active,                            &
+                  &               deeppos_active, nutprof_active,                            &
+                  &               bstate_2d, outincs_2d,                                     &
+                  &               diag_active, diag,                                         &
+                  &               diag_fulldepth_active, diag_fulldepth_2d )
+
+               outincs_3d(:,:,jk,:) = outincs_2d(:,:,1,:)
+            END DO
+         ENDIF
+
+         ! Save balancing increments
+         phyto_balinc(:,:,:,jp_had_nut) = outincs_3d(:,:,:,i_tracer(1))
+         phyto_balinc(:,:,:,jp_had_phy) = outincs_3d(:,:,:,i_tracer(2))
+         phyto_balinc(:,:,:,jp_had_zoo) = outincs_3d(:,:,:,i_tracer(3))
+         phyto_balinc(:,:,:,jp_had_pdn) = outincs_3d(:,:,:,i_tracer(4))
+         phyto_balinc(:,:,:,jp_had_dic) = outincs_3d(:,:,:,i_tracer(5))
+         phyto_balinc(:,:,:,jp_had_alk) = outincs_3d(:,:,:,i_tracer(6))
+      
+      ELSE   ! No nitrogen balancing
+      
+         ! Initialise phytoplankton increment to zero
+         phyto_balinc(:,:,:,jp_had_phy) = 0.0
+         
+         ! Convert surface chlorophyll increment to phytoplankton nitrogen
+         DO jk = 1, kdeps
+            phyto_balinc(:,:,jk,jp_had_phy) = ( cchl_p_bkg_3d(:,:,jk) / (mw_carbon * c2n_p) ) * pinc_chltot_3d(:,:,jk)
+         END DO
+         
+         IF (kdeps == 1) THEN
+            ! Propagate through mixed layer
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !
+                  jkmax = jpk-1
+                  DO jk = jpk-1, 1, -1
+                     IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                        & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                        pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                        jkmax = jk
+                     ENDIF
+                  END DO
+                  !
+                  DO jk = 2, jkmax
+                     phyto_balinc(ji,jj,jk,jp_had_phy) = phyto_balinc(ji,jj,1,jp_had_phy)
+                  END DO
+                  !
+               END DO
+            END DO
+         ENDIF
+
+         ! Set other balancing increments to zero
+         phyto_balinc(:,:,:,jp_had_nut) = 0.0
+         phyto_balinc(:,:,:,jp_had_zoo) = 0.0
+         phyto_balinc(:,:,:,jp_had_pdn) = 0.0
+         phyto_balinc(:,:,:,jp_had_dic) = 0.0
+         phyto_balinc(:,:,:,jp_had_alk) = 0.0
+      
+      ENDIF
+
+   END SUBROUTINE asm_phyto_bal_hadocc
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option : Empty routine
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE asm_phyto_bal_hadocc( kdeps,                              &
+      &                             ld_chltot,                          &
+      &                             pinc_chltot_3d,                     &
+      &                             ld_phytot,                          &
+      &                             pinc_phytot,                        &
+      &                             pincper,                            &
+      &                             p_maxchlinc, ld_phytobal, pmld,     &
+      &                             pgrow_avg_bkg_3d, ploss_avg_bkg_3d, &
+      &                             phyt_avg_bkg_3d, mld_max_bkg,       &
+      &                             cchl_p_bkg_3d,                      &
+      &                             tracer_bkg, phyto_balinc )
+      INTEGER :: kdeps
+      LOGICAL :: ld_chltot
+      REAL    :: pinc_chltot_3d(:,:,:)
+      LOGICAL :: ld_phytot
+      REAL    :: pinc_phytot_3d(:,:,:)
+      REAL    :: pincper
+      REAL    :: p_maxchlinc
+      LOGICAL :: ld_phytobal
+      REAL    :: pmld(:,:)
+      REAL    :: pgrow_avg_bkg_3d(:,:,:)
+      REAL    :: ploss_avg_bkg_3d(:,:,:)
+      REAL    :: phyt_avg_bkg_3d(:,:,:)
+      REAL    :: mld_max_bkg(:,:)
+      REAL    :: cchl_p_bkg_3d(:,:,:)
+      REAL    :: tracer_bkg(:,:,:,:)
+      REAL    :: phyto_balinc(:,:,:,:)
+      WRITE(*,*) 'asm_phyto_bal_hadocc: You should not have seen this print! error?'
+   END SUBROUTINE asm_phyto_bal_hadocc
+#endif
+
+   !!======================================================================
+END MODULE asmphytobal_hadocc
Index: /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_medusa.F90
===================================================================
--- /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_medusa.F90	(revision 13354)
+++ /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_medusa.F90	(revision 13354)
@@ -0,0 +1,568 @@
+MODULE asmphytobal_medusa
+   !!======================================================================
+   !!                       ***  MODULE asmphyto2dbal_medusa  ***
+   !! Calculate increments to MEDUSA based on surface phyto2d increments
+   !!
+   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
+   !! For licensing reasons this is kept in its own internal Met Office
+   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
+   !! repository, and must be merged in when building.
+   !!
+   !!======================================================================
+   !! History :  3.6  ! 2017-08 (D. Ford)  Adapted from asmphyto2dbal_hadocc
+   !!----------------------------------------------------------------------
+#if defined key_asminc && defined key_medusa
+   !!----------------------------------------------------------------------
+   !! 'key_asminc'          : assimilation increment interface
+   !! 'key_medusa'          : MEDUSA model
+   !!----------------------------------------------------------------------
+   !! asm_phyto2d_bal_medusa : routine to calculate increments to MEDUSA
+   !!----------------------------------------------------------------------
+   USE par_kind,      ONLY: wp             ! kind parameters
+   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
+   USE dom_oce,       ONLY: gdepw_n        ! domain information
+   USE zdftmx,        ONLY: ln_tmx_itf, &  ! Indonesian Throughflow
+      &                     mask_itf       ! tidal mixing mask
+   USE iom                                 ! i/o
+   USE sms_medusa                          ! MEDUSA parameters
+   USE par_medusa                          ! MEDUSA parameters
+   USE par_trc,       ONLY: jptra          ! Tracer parameters
+   USE bioanalysis                         ! Nitrogen balancing
+
+   IMPLICIT NONE
+   PRIVATE                   
+
+   PUBLIC asm_phyto_bal_medusa
+
+   ! Default values for biological assimilation parameters
+   ! Should match Hemmings et al. (2008)
+   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
+   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
+   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
+   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
+   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
+   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
+   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
+   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
+   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
+   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
+   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
+   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
+   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
+   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
+   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
+   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
+   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
+   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
+   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
+   !
+   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
+   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
+   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
+   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
+   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
+   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
+   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
+   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
+
+CONTAINS
+
+   SUBROUTINE asm_phyto_bal_medusa( kdeps,                              &
+      &                             ld_chltot,                          &
+      &                             pinc_chltot_3d,                     &
+      &                             ld_chldia,                          &
+      &                             pinc_chldia_3d,                     &
+      &                             ld_chlnon,                          &
+      &                             pinc_chlnon_3d,                     &
+      &                             ld_phytot,                          &
+      &                             pinc_phytot_3d,                     &
+      &                             ld_phydia,                          &
+      &                             pinc_phydia_3d,                     &
+      &                             ld_phynon,                          &
+      &                             pinc_phynon_3d,                     &
+      &                             pincper,                            &
+      &                             p_maxchlinc, ld_phytobal, pmld,     &
+      &                             pgrow_avg_bkg_3d, ploss_avg_bkg_3d, &
+      &                             phyt_avg_bkg_3d, mld_max_bkg,       &
+      &                             tracer_bkg, phyto_balinc )
+      !!---------------------------------------------------------------------------
+      !!                    ***  ROUTINE asm_phyto_bal_medusa  ***
+      !!
+      !! ** Purpose :   calculate increments to MEDUSA from phytoplankton increments
+      !!
+      !! ** Method  :   average up MEDUSA to look like HadOCC
+      !!                call nitrogen balancing scheme
+      !!                separate back out to MEDUSA
+      !!
+      !! ** Action  :   populate phyto_balinc
+      !!
+      !! References :   Hemmings et al., 2008, J. Mar. Res.
+      !!                Ford et al., 2012, Ocean Sci.
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER,  INTENT(in   )                               :: kdeps            ! No. inc deps 1 or jpk
+      LOGICAL,  INTENT(in   )                               :: ld_chltot        ! Assim chltot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_chltot_3d   ! chltot increments (3D)
+      LOGICAL,  INTENT(in   )                               :: ld_chldia        ! Assim chldia y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_chldia_3d   ! chldia increments (3D)
+      LOGICAL,  INTENT(in   )                               :: ld_chlnon        ! Assim chlnon y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_chlnon_3d   ! chlnon increments (3D)
+      LOGICAL,  INTENT(in   )                               :: ld_phytot        ! Assim phytot y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_phytot_3d   ! phytot increments (3D)
+      LOGICAL,  INTENT(in   )                               :: ld_phydia        ! Assim phydia y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_phydia_3d   ! phydia increments (3D)
+      LOGICAL,  INTENT(in   )                               :: ld_phynon        ! Assim phynon y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_phynon_3d   ! phynon increments (3D)
+      REAL(wp), INTENT(in   )                               :: pincper          ! Assimilation period
+      REAL(wp), INTENT(in   )                               :: p_maxchlinc      ! Max chl increment
+      LOGICAL,  INTENT(in   )                               :: ld_phytobal      ! Balancing y/n
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld             ! Mixed layer depth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: pgrow_avg_bkg_3d ! Avg phyto growth (3D)
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: ploss_avg_bkg_3d ! Avg phyto loss (3D)
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: phyt_avg_bkg_3d  ! Avg phyto (3D)
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg      ! Max MLD
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg       ! State variables
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto_balinc     ! Balancing increments
+      !!
+      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
+      INTEGER                                               :: jkmax          ! Loop index
+      INTEGER                                               :: jkinc          ! Loop index
+      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
+      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
+      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
+      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
+      REAL(wp)                                              :: zfrac          ! 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),                DIMENSION(jpi,jpj)           :: cchl_p_2d      ! C:Chl for total phy (2D)
+      REAL(wp),                DIMENSION(jpi,jpj,jpk)       :: cchl_p_3d      ! C:Chl for total phy (3D)
+      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
+      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
+      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: bstate_2d      ! Background state (2D)
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate_3d      ! Background state (3D)
+      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: outincs_2d     ! Balancing increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs_3d     ! Balancing increments (3D)
+      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag              ! Depth-indep diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj,1,22)      :: diag_fulldepth_2d ! Full-depth diagnostics (2D)
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth_3d ! Full-depth diagnostics (3D)
+      REAL(wp),                DIMENSION(jpi,jpj,1)         :: tmask_2d          ! Single-level tmask
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_2d    ! chltot increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chldia_2d    ! chldia increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chlnon_2d    ! chlnon increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_phytot_2d    ! phytot increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_phydia_2d    ! phydia increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_phynon_2d    ! phynon increments (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pgrow_avg_bkg_2d  ! Avg phyto growth (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: ploss_avg_bkg_2d  ! Avg phyto loss (2D)
+      REAL(wp),                DIMENSION(jpi,jpj)           :: phyt_avg_bkg_2d   ! Avg phyto (2D)
+      !!---------------------------------------------------------------------------
+
+      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      IF ( p_maxchlinc > 0.0 ) THEN
+         DO jk = 1, kdeps
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( ld_chltot ) THEN
+                     pinc_chltot_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_3d(ji,jj,jk), p_maxchlinc ) )
+                  ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
+                     pinc_chltot_3d(ji,jj,jk) = pinc_chldia_3d(ji,jj,jk) + pinc_chlnon_3d(ji,jj,jk)
+                     pinc_chltot_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_3d(ji,jj,jk), p_maxchlinc ) )
+                     IF ( pinc_chltot_3d(ji,jj,jk) .NE. ( pinc_chldia_3d(ji,jj,jk) + pinc_chlnon_3d(ji,jj,jk) ) ) THEN
+                        zfrac = pinc_chltot_3d(ji,jj,jk) / ( pinc_chldia_3d(ji,jj,jk) + pinc_chlnon_3d(ji,jj,jk) )
+                        pinc_chldia_3d(ji,jj,jk) = pinc_chldia_3d(ji,jj,jk) * zfrac
+                        pinc_chlnon_3d(ji,jj,jk) = pinc_chlnon_3d(ji,jj,jk) * zfrac
+                     ENDIF
+                  ELSE IF ( ld_chldia ) THEN
+                     pinc_chldia_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chldia_3d(ji,jj,jk), p_maxchlinc ) )
+                     pinc_chltot_3d(ji,jj,jk) = pinc_chldia_3d(ji,jj,jk)
+                  ELSE IF ( ld_chlnon ) THEN
+                     pinc_chlnon_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chlnon_3d(ji,jj,jk), p_maxchlinc ) )
+                     pinc_chltot_3d(ji,jj,jk) = pinc_chlnon_3d(ji,jj,jk)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      IF ( ld_phytot .OR. ld_phydia .OR. ld_phynon ) THEN
+         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
+      ENDIF
+      
+      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
+
+         ! Set up model parameters to be passed into Hemmings balancing routine.
+         ! For now these are hardwired to the standard HadOCC parameter values
+         ! (except C:N ratios) as this is what the scheme was developed for.
+         ! Obviously, HadOCC and MEDUSA are rather different models, so this
+         ! isn't ideal, but there's not always direct analogues between the two
+         ! parameter sets, so it's the easiest way to get something running.
+         ! In the longer term, some serious MarMOT-based development is required.
+         modparm(1)  = 0.1                               ! grow_sat
+         modparm(2)  = 2.0                               ! psmax
+         modparm(3)  = 0.845                             ! par
+         modparm(4)  = 0.02                              ! alpha
+         modparm(5)  = 0.05                              ! resp_rate
+         modparm(6)  = 0.05                              ! pmort_rate
+         modparm(7)  = 0.01                              ! phyto_min
+         modparm(8)  = 0.05                              ! z_mort_1
+         modparm(9)  = 1.0                               ! z_mort_2
+         modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p
+         modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z
+         modparm(12) = xthetad                           ! c2n_d
+         modparm(13) = 0.01                              ! graze_threshold
+         modparm(14) = 2.0                               ! holling_coef
+         modparm(15) = 0.5                               ! graze_sat
+         modparm(16) = 2.0                               ! graze_max
+
+         ! Set up assimilation parameters to be passed into balancing routine
+         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
+         assimparm(2)  = balnutext
+         assimparm(3)  = balnutmin
+         assimparm(4)  = r
+         assimparm(5)  = beta_g
+         assimparm(6)  = beta_l
+         assimparm(7)  = beta_m
+         assimparm(8)  = a_g
+         assimparm(9)  = a_l
+         assimparm(10) = a_m
+         assimparm(11) = zfracb0
+         assimparm(12) = zfracb1
+         assimparm(13) = qrfmax
+         assimparm(14) = qafmax
+         assimparm(15) = zrfmax
+         assimparm(16) = zafmax
+         assimparm(17) = prfmax
+         assimparm(18) = incphymin
+         assimparm(19) = integnstep
+         assimparm(20) = pthreshold
+
+         ! Set up external tracer indices array bstate
+         i_tracer(1) = 1   ! nutrient
+         i_tracer(2) = 2   ! phytoplankton
+         i_tracer(3) = 3   ! zooplankton
+         i_tracer(4) = 4   ! detritus
+         i_tracer(5) = 5   ! DIC
+         i_tracer(6) = 6   ! Alkalinity
+
+         ! Set background state
+         bstate_3d(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin)
+         bstate_3d(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd)
+         bstate_3d(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme)
+         bstate_3d(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet)
+         bstate_3d(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic)
+         bstate_3d(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk)
+
+         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
+         ! and nitrogen to biomass equivalent for PZD
+         ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA
+         cchl_p_3d(:,:,:) = 0.0
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF ( ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd ) ) .GT. 0.0 ) THEN
+                     cchl_p_3d(ji,jj,jk) = xmassc * ( ( tracer_bkg(ji,jj,jk,jpphn) * xthetapn ) +      &
+                        &                             ( tracer_bkg(ji,jj,jk,jpphd) * xthetapd )   ) /  &
+                        &                  ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd ) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) )
+         n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) )
+         n2be_d = 14.01 + ( xmassc * xthetad )
+
+         ! Call nitrogen balancing routine
+         IF (kdeps == 1) THEN
+            pinc_chltot_2d(:,:)   = pinc_chltot_3d(:,:,1)
+            cchl_p_2d(:,:)        = cchl_p_3d(:,:,1)
+            phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg_3d(:,:,1)
+            pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg_3d(:,:,1)
+            ploss_avg_bkg_2d(:,:) = ploss_avg_bkg_3d(:,:,1)
+            
+            CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
+               &               n2be_p, n2be_z, n2be_d, assimparm,                      &
+               &               INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
+               &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot_2d(:,:), cchl_p_2d(:,:), &
+               &               nbal_active, phyt_avg_bkg_2d(:,:),                      &
+               &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:), &
+               &               subsurf_active, deepneg_active,                         &
+               &               deeppos_active, nutprof_active,                         &
+               &               bstate_3d, outincs_3d,                                  &
+               &               diag_active, diag,                                      &
+               &               diag_fulldepth_active, diag_fulldepth_3d )
+         ELSE
+            pmld(:,:) = 0.5
+            
+            DO jk = 1, kdeps
+               pinc_chltot_2d(:,:)   = pinc_chltot_3d(:,:,jk)
+               cchl_p_2d(:,:)        = cchl_p_3d(:,:,jk)
+               phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg_3d(:,:,jk)
+               pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg_3d(:,:,jk)
+               ploss_avg_bkg_2d(:,:) = ploss_avg_bkg_3d(:,:,jk)
+               tmask_2d(:,:,1)       = tmask(:,:,jk)
+               bstate_2d(:,:,1,:)    = bstate_3d(:,:,jk,:)
+               outincs_2d(:,:,:,:)   = 0.0
+
+               CALL bio_analysis( jpi, jpj, 1, gdepw_n(:,:,2), i_tracer, modparm,            &
+                  &               n2be_p, n2be_z, n2be_d, assimparm,                         &
+                  &               INT(pincper), 1, INT(SUM(tmask_2d,3)), tmask_2d(:,:,:),    &
+                  &               pmld(:,:), pmld(:,:), pinc_chltot_2d(:,:), cchl_p_2d(:,:), &
+                  &               nbal_active, phyt_avg_bkg_2d(:,:),                         &
+                  &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:),   &
+                  &               subsurf_active, deepneg_active,                            &
+                  &               deeppos_active, nutprof_active,                            &
+                  &               bstate_2d, outincs_2d,                                     &
+                  &               diag_active, diag,                                         &
+                  &               diag_fulldepth_active, diag_fulldepth_2d )
+
+               outincs_3d(:,:,jk,:) = outincs_2d(:,:,1,:)
+            END DO
+         ENDIF
+         
+         ! Loop over each grid point partioning the increments
+         phyto_balinc(:,:,:,:) = 0.0
+         DO jk = 1, jpk
+            IF (kdeps == 1) THEN
+               jkinc = 1
+            ELSE
+               IF (jk > kdeps) THEN
+                  EXIT
+               ENDIF
+               jkinc = jk
+            ENDIF
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  ! Phytoplankton
+                  IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. &
+                     & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. &
+                     & ( pinc_chltot_3d(ji,jj,jkinc) /= 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
+                        ! Phytoplankton nitrogen split up based on assimilation increments
+                        zfrac_phn = pinc_chlnon_3d(ji,jj,jkinc) / pinc_chltot_3d(ji,jj,jkinc)
+                        zfrac_phd = pinc_chldia_3d(ji,jj,jkinc) / pinc_chltot_3d(ji,jj,jkinc)
+                     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)
+                     
+                     phyto_balinc(ji,jj,jk,jpphn) = outincs_3d(ji,jj,jk,i_tracer(2)) * zfrac_phn
+                     phyto_balinc(ji,jj,jk,jpphd) = outincs_3d(ji,jj,jk,i_tracer(2)) * zfrac_phd
+                     phyto_balinc(ji,jj,jk,jppds) = phyto_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
+                     phyto_balinc(ji,jj,jk,jpchn) = phyto_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
+                     phyto_balinc(ji,jj,jk,jpchd) = phyto_balinc(ji,jj,jk,jpphd) * zrat_chd_phd
+                  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))
+                     phyto_balinc(ji,jj,jk,jpzmi) = outincs_3d(ji,jj,jk,i_tracer(3)) * zfrac_zmi
+                     phyto_balinc(ji,jj,jk,jpzme) = outincs_3d(ji,jj,jk,i_tracer(3)) * zfrac_zme
+                  ENDIF
+
+                  ! Nitrogen nutrient straight from balancing scheme
+                  phyto_balinc(ji,jj,jk,jpdin) = outincs_3d(ji,jj,jk,i_tracer(1))
+
+                  ! Nitrogen detritus straight from balancing scheme
+                  phyto_balinc(ji,jj,jk,jpdet) = outincs_3d(ji,jj,jk,i_tracer(4))
+
+                  ! DIC straight from balancing scheme
+                  phyto_balinc(ji,jj,jk,jpdic) = outincs_3d(ji,jj,jk,i_tracer(5))
+
+                  ! Alkalinity straight from balancing scheme
+                  phyto_balinc(ji,jj,jk,jpalk) = outincs_3d(ji,jj,jk,i_tracer(6))
+
+                  ! Remove diatom silicate increment from nutrient silicate to conserve mass
+                  IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
+                     phyto_balinc(ji,jj,jk,jpsil) = phyto_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)
+                     phyto_balinc(ji,jj,jk,jpdtc) = phyto_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
+                  ENDIF
+
+                  ! Do nothing with iron or oxygen for the time being
+                  phyto_balinc(ji,jj,jk,jpfer) = 0.0
+                  phyto_balinc(ji,jj,jk,jpoxy) = 0.0
+                  
+               END DO
+            END DO
+         END DO
+      
+      ELSE   ! No nitrogen balancing
+      
+         ! Initialise individual chlorophyll increments to zero
+         phyto_balinc(:,:,:,jpchn) = 0.0
+         phyto_balinc(:,:,:,jpchd) = 0.0
+         
+         ! Split up total surface chlorophyll increments
+         DO jk = 1, kdeps
+            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
+                     IF ( ld_chltot ) THEN
+                        ! Chlorophyll split up based on existing ratios
+                        zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / &
+                           &        ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd) )
+                        zfrac_chd = tracer_bkg(ji,jj,jk,jpchd) / &
+                           &        ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd) )
+                        phyto_balinc(ji,jj,jk,jpchn) = pinc_chltot_3d(ji,jj,jk) * zfrac_chn
+                        phyto_balinc(ji,jj,jk,jpchd) = pinc_chltot_3d(ji,jj,jk) * zfrac_chd
+                     ENDIF
+                     IF( ld_chldia ) THEN
+                        phyto_balinc(ji,jj,jk,jpchd) = pinc_chldia_3d(ji,jj,jk)
+                     ENDIF
+                     IF( ld_chlnon ) THEN
+                        phyto_balinc(ji,jj,jk,jpchn) = pinc_chlnon_3d(ji,jj,jk)
+                     ENDIF
+
+                     ! Maintain stoichiometric ratios of nitrogen and silicate
+                     IF ( ld_chltot .OR. ld_chlnon ) THEN
+                        zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
+                        phyto_balinc(ji,jj,jk,jpphn) = phyto_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
+                     ENDIF
+                     IF ( ld_chltot .OR. ld_chldia ) THEN
+                        zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
+                        phyto_balinc(ji,jj,jk,jpphd) = phyto_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
+                        zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
+                        phyto_balinc(ji,jj,jk,jppds) = phyto_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
+                     ENDIF
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         
+         IF (kdeps == 1) THEN
+            ! Propagate through mixed layer
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !
+                  jkmax = jpk-1
+                  DO jk = jpk-1, 1, -1
+                     IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
+                        & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
+                        pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
+                        jkmax = jk
+                     ENDIF
+                  END DO
+                  !
+                  DO jk = 2, jkmax
+                     phyto_balinc(ji,jj,jk,jpchn) = phyto_balinc(ji,jj,1,jpchn)
+                     phyto_balinc(ji,jj,jk,jpchd) = phyto_balinc(ji,jj,1,jpchd)
+                     phyto_balinc(ji,jj,jk,jpphn) = phyto_balinc(ji,jj,1,jpphn)
+                     phyto_balinc(ji,jj,jk,jpphd) = phyto_balinc(ji,jj,1,jpphd)
+                     phyto_balinc(ji,jj,jk,jppds) = phyto_balinc(ji,jj,1,jppds)
+                  END DO
+                  !
+               END DO
+            END DO
+         ENDIF
+
+         ! Set other balancing increments to zero
+         phyto_balinc(:,:,:,jpzmi) = 0.0
+         phyto_balinc(:,:,:,jpzme) = 0.0
+         phyto_balinc(:,:,:,jpdin) = 0.0
+         phyto_balinc(:,:,:,jpsil) = 0.0
+         phyto_balinc(:,:,:,jpfer) = 0.0
+         phyto_balinc(:,:,:,jpdet) = 0.0
+         phyto_balinc(:,:,:,jpdtc) = 0.0
+         phyto_balinc(:,:,:,jpdic) = 0.0
+         phyto_balinc(:,:,:,jpalk) = 0.0
+         phyto_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
+               phyto_balinc(:,:,jk,jn) = phyto_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
+            END DO
+         END DO
+      ENDIF
+
+   END SUBROUTINE asm_phyto_bal_medusa
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option : Empty routine
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE asm_phyto_bal_medusa( kdeps,                          &
+      &                             ld_chltot,                      &
+      &                             pinc_chltot_3d,                    &
+      &                             ld_chldia,                      &
+      &                             pinc_chldia_3d,                    &
+      &                             ld_chlnon,                      &
+      &                             pinc_chlnon_3d,                    &
+      &                             ld_phytot,                      &
+      &                             pinc_phytot_3d,                    &
+      &                             ld_phydia,                      &
+      &                             pinc_phydia_3d,                    &
+      &                             ld_phynon,                      &
+      &                             pinc_phynon_3d,                    &
+      &                             pincper,                        &
+      &                             p_maxchlinc, ld_phytobal, pmld, &
+      &                             pgrow_avg_bkg_3d, ploss_avg_bkg_3d,   &
+      &                             phyt_avg_bkg_3d, mld_max_bkg,      &
+      &                             tracer_bkg, phyto_balinc )
+      INTEGER :: kdeps
+      LOGICAL :: ld_chltot
+      REAL    :: pinc_chltot_3d(:,:,:)
+      LOGICAL :: ld_chldia
+      REAL    :: pinc_chldia_3d(:,:,:)
+      LOGICAL :: ld_chlnon
+      REAL    :: pinc_chlnon_3d(:,:,:)
+      LOGICAL :: ld_phytot
+      REAL    :: pinc_phytot_3d(:,:,:)
+      LOGICAL :: ld_phydia
+      REAL    :: pinc_phydia_3d(:,:,:)
+      LOGICAL :: ld_phynon
+      REAL    :: pinc_phynon_3d(:,:,:)
+      REAL    :: pincper
+      REAL    :: p_maxchlinc
+      LOGICAL :: ld_phytobal
+      REAL    :: pmld(:,:)
+      REAL    :: pgrow_avg_bkg_3d(:,:,:)
+      REAL    :: ploss_avg_bkg_3d(:,:,:)
+      REAL    :: phyt_avg_bkg_3d(:,:,:)
+      REAL    :: mld_max_bkg(:,:)
+      REAL    :: tracer_bkg(:,:,:,:)
+      REAL    :: phyto_balinc(:,:,:,:)
+      WRITE(*,*) 'asm_phyto_bal_medusa: You should not have seen this print! error?'
+   END SUBROUTINE asm_phyto_bal_medusa
+#endif
+
+   !!======================================================================
+END MODULE asmphytobal_medusa
Index: /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 13353)
+++ /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 13354)
@@ -35,5 +35,6 @@
       USE bio_medusa_mod
       USE par_oce,           ONLY: jpi, jpj, jpk
-      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max
+      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max, &
+         &                         pgrow_avg_3d, ploss_avg_3d, phyt_avg_3d
       USE trc,               ONLY: ln_diatrc, med_diag, nittrc000 
       USE in_out_manager,    ONLY: lwp, numout
@@ -199,4 +200,7 @@
       ploss_avg(:,:) = 0.0
       phyt_avg(:,:)  = 0.0
+      pgrow_avg_3d(:,:,:) = 0.0
+      ploss_avg_3d(:,:,:) = 0.0
+      phyt_avg_3d(:,:,:)  = 0.0
       IF( kt == nittrc000 ) THEN
          mld_max(:,:) = 0.0
Index: /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 13353)
+++ /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 13354)
@@ -46,4 +46,5 @@
                                    ln_foam_medusa,                         &
                                    pgrow_avg, ploss_avg, phyt_avg,         &
+                                   pgrow_avg_3d, ploss_avg_3d, phyt_avg_3d, &
                                    xkphd, xkphn, xkzme, xkzmi,             &
                                    xmetapd, xmetapn, xmetazme, xmetazmi,   &
@@ -229,4 +230,12 @@
                                      ((zphn(ji,jj) + zphd(ji,jj)) *       &
                                       fse3t(ji,jj,jk) * fq0)
+                  !!
+                  pgrow_avg_3d(ji,jj,jk) = (fprn(ji,jj) * zphn(ji,jj)) +  &
+                                           (fprd(ji,jj) * zphd(ji,jj))
+                  ploss_avg_3d(ji,jj,jk) = fgmepd(ji,jj) + fdpd(ji,jj) +  &
+                                           fdpd2(ji,jj)                +  &
+                                           fgmepn(ji,jj) + fdpn(ji,jj) +  &
+                                           fdpn2(ji,jj)  + fgmipn(ji,jj)
+                  phyt_avg_3d(ji,jj,jk)  = zphn(ji,jj) + zphd(ji,jj)
                ENDIF
             ENDDO
Index: /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90
===================================================================
--- /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 13353)
+++ /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 13354)
@@ -362,4 +362,7 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ploss_avg  !: Mixed layer average phytoplankton loss
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: phyt_avg   !: Mixed layer average phytoplankton
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: pgrow_avg_3d  !: Mixed layer average phytoplankton growth
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ploss_avg_3d  !: Mixed layer average phytoplankton loss
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: phyt_avg_3d   !: Mixed layer average phytoplankton
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_max    !: Maximum mixed layer depth
 !!
@@ -438,4 +441,6 @@
       !* Fields for ocean colour data assimilation
       ALLOCATE( pgrow_avg(jpi,jpj)   , ploss_avg(jpi,jpj)   ,       &
+         &      pgrow_avg_3d(jpi,jpj,jpk) , ploss_avg_3d(jpi,jpj,jpk) , &
+         &      phyt_avg_3d(jpi,jpj,jpk)  ,                         &
          &      phyt_avg(jpi,jpj)    , mld_max(jpi,jpj)     ,    STAT=ierr(9) )
 #endif
Index: /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 13353)
+++ /branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 13354)
@@ -368,4 +368,15 @@
             mld_max(:,:)   = 0.0
          ENDIF
+         IF( iom_varid( numrtr, 'pgrow_avg_3d', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg present - reading in ...'
+            CALL iom_get( numrtr, jpdom_autoglo, 'pgrow_avg_3d',  pgrow_avg_3d(:,:,:)  )
+            CALL iom_get( numrtr, jpdom_autoglo, 'ploss_avg_3d',  ploss_avg_3d(:,:,:)  )
+            CALL iom_get( numrtr, jpdom_autoglo, 'phyt_avg_3d',   phyt_avg_3d(:,:,:)   )
+         ELSE
+            IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg_3d absent - setting to zero ...'
+            pgrow_avg_3d(:,:,:) = 0.0
+            ploss_avg_3d(:,:,:) = 0.0
+            phyt_avg_3d(:,:,:)  = 0.0
+         ENDIF
       ENDIF
 
@@ -553,4 +564,7 @@
          CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg',  phyt_avg(:,:)  )
          CALL iom_rstput( kt, nitrst, numrtw, 'mld_max',   mld_max(:,:)   )
+         CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg_3d', pgrow_avg_3d(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg_3d', ploss_avg_3d(:,:,:) )
+         CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg_3d',  phyt_avg_3d(:,:,:)  )
       ENDIF
 !!
