Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 11989)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90	(revision 11990)
@@ -59,4 +59,6 @@
    USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
       & asm_phyto2d_bal_medusa
+   USE asmphyto3dbal_medusa, ONLY: & ! phyto3d balancing for MEDUSA
+      & asm_phyto3d_bal_medusa
    USE asmpco2bal, ONLY:    & ! pCO2 balancing for MEDUSA
       & asm_pco2_bal
@@ -71,4 +73,7 @@
       & ploss_avg,          &
       & phyt_avg,           &
+      & pgrow_avg_3d,       &
+      & ploss_avg_3d,       &
+      & phyt_avg_3d,        &
       & mld_max
 #elif defined key_hadocc
@@ -172,4 +177,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
@@ -213,7 +221,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.
@@ -525,4 +534,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) )
@@ -530,4 +542,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
@@ -565,8 +580,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)
 
@@ -727,4 +748,16 @@
             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) )
@@ -810,4 +843,7 @@
       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          )
       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) )
@@ -1163,4 +1199,5 @@
       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
@@ -1215,25 +1252,35 @@
             END DO
          ENDIF
-
+         
 #if defined key_medusa && defined key_foam_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
+         IF ( ln_phytobal ) THEN
+            zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
+            CALL asm_phyto3d_bal_medusa( chl_inc,                             &
+               &                         zincper,                             &
+               &                         rn_maxchlinc,                        &
+               &                         pgrow_avg_3d_bkg, ploss_avg_3d_bkg,  &
+               &                         phyt_avg_3d_bkg,                     &
+               &                         tracer_bkg, phyto3d_balinc )
+         ELSE
+            ! 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
-         END DO
+         ENDIF
 #elif defined key_hadocc
          phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
@@ -1420,6 +1467,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 ( ln_slchltotinc ) 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 ( ln_plchltotinc ) 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)
@@ -1666,8 +1719,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 ( ln_slchltotinc ) 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 ( ln_plchltotinc ) 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)
@@ -1861,5 +1922,5 @@
             CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
 #endif
-            IF ( ln_phytobal ) THEN
+            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN
                pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
@@ -1881,5 +1942,5 @@
             CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
 #endif
-            IF ( ln_phytobal ) THEN
+            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN
                psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
@@ -1903,5 +1964,5 @@
             CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
 #endif
-            IF ( ln_phytobal ) THEN
+            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN
                pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
@@ -1925,5 +1986,5 @@
             CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
 #endif
-            IF ( ln_phytobal ) THEN
+            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN
                palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
@@ -1945,5 +2006,5 @@
             CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
 #endif
-            IF ( ln_phytobal ) THEN
+            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN
                po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
             ENDIF
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto3dbal_medusa.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto3dbal_medusa.F90	(revision 11990)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto3dbal_medusa.F90	(revision 11990)
@@ -0,0 +1,373 @@
+MODULE asmphyto3dbal_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 && defined key_foam_medusa
+   !!----------------------------------------------------------------------
+   !! 'key_asminc'          : assimilation increment interface
+   !! 'key_medusa'          : MEDUSA model
+   !! 'key_foam_medusa'     : MEDUSA extras for FOAM OBS and ASM
+   !!----------------------------------------------------------------------
+   !! 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_phyto3d_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_phyto3d_bal_medusa( pinc_chltot,                    &
+      &                               pincper,                        &
+      &                               p_maxchlinc,                    &
+      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                               phyt_avg_bkg,                   &
+      &                               tracer_bkg, phyto3d_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.
+      !!---------------------------------------------------------------------------
+      !!
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk)       :: pinc_chltot    ! chltot increments
+      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
+      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: pgrow_avg_bkg  ! Avg phyto growth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: ploss_avg_bkg  ! Avg phyto loss
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: phyt_avg_bkg   ! Avg phyto
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto3d_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,jpk)       :: cchl_p         ! C:Chl for total phy
+      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p_2d      ! 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,1,6)       :: bstate_2d      ! Background state
+      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
+      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: outincs_2d     ! Balancing increments
+      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj,1,22)      :: diag_fulldepth ! Full-depth diagnostics
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_2d
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pgrow_avg_bkg_2d
+      REAL(wp),                DIMENSION(jpi,jpj)           :: ploss_avg_bkg_2d
+      REAL(wp),                DIMENSION(jpi,jpj)           :: phyt_avg_bkg_2d
+      REAL(wp),                DIMENSION(jpi,jpj,1)         :: tmask_2d
+      REAL(wp),                DIMENSION(jpi,jpj)           :: pmld_2d
+      REAL(wp),                DIMENSION(jpi,jpj)           :: mld_max_bkg_2d
+      !!---------------------------------------------------------------------------
+
+      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      IF ( p_maxchlinc > 0.0 ) THEN
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  pinc_chltot(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj,jk), p_maxchlinc ) )
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      ! 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 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(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 )
+
+      pmld_2d(:,:)        = 0.5
+      mld_max_bkg_2d(:,:) = 0.5
+
+      DO jk = 1, jpk
+
+         tmask_2d(:,:,1)       = tmask(:,:,jk)
+         pinc_chltot_2d(:,:)   = pinc_chltot(:,:,jk)
+         cchl_p_2d(:,:)        = cchl_p(:,:,jk)
+         phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg(:,:,jk)
+         pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg(:,:,jk)
+         ploss_avg_bkg_2d(:,:) = ploss_avg_bkg(:,:,jk)
+         bstate_2d(:,:,1,:)    = bstate(:,:,jk,:)
+         outincs_2d(:,:,:,:)   = 0.0
+         
+         ! Call nitrogen balancing routine
+         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_2d(:,:), mld_max_bkg_2d(:,:), 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 )
+
+         outincs(:,:,jk,:) = outincs_2d(:,:,1,:)
+
+      END DO
+
+      ! Loop over each grid point partioning the increments
+      phyto3d_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,jk) /= 0.0 ) ) 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))
+
+                  ! 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)
+
+                  phyto3d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn
+                  phyto3d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd
+                  phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
+                  phyto3d_balinc(ji,jj,jk,jpchn) = phyto3d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
+                  phyto3d_balinc(ji,jj,jk,jpchd) = phyto3d_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))
+                  phyto3d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi
+                  phyto3d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme
+               ENDIF
+
+               ! Nitrogen nutrient straight from balancing scheme
+               phyto3d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1))
+
+               ! Nitrogen detritus straight from balancing scheme
+               phyto3d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4))
+
+               ! DIC straight from balancing scheme
+               phyto3d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5))
+
+               ! Alkalinity straight from balancing scheme
+               phyto3d_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) - phyto3d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
+                  phyto3d_balinc(ji,jj,jk,jpsil) = phyto3d_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)
+                  phyto3d_balinc(ji,jj,jk,jpdtc) = phyto3d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
+               ENDIF
+
+               ! Do nothing with iron or oxygen for the time being
+               phyto3d_balinc(ji,jj,jk,jpfer) = 0.0
+               phyto3d_balinc(ji,jj,jk,jpoxy) = 0.0
+
+            END DO
+         END DO
+      END DO
+      
+      ! 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
+               phyto3d_balinc(:,:,jk,jn) = phyto3d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
+            END DO
+         END DO
+      ENDIF
+
+   END SUBROUTINE asm_phyto3d_bal_medusa
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option : Empty routine
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE asm_phyto3d_bal_medusa(pinc_chltot,                    &
+      &                              pincper,                        &
+      &                              p_maxchlinc,                    &
+      &                              pgrow_avg_bkg, ploss_avg_bkg,   &
+      &                              phyt_avg_bkg,                   &
+      &                              tracer_bkg, phyto3d_balinc )
+      REAL    :: pinc_chltot(:,:,:)
+      REAL    :: pincper
+      REAL    :: p_maxchlinc
+      REAL    :: pgrow_avg_bkg(:,:,:)
+      REAL    :: ploss_avg_bkg(:,:,:)
+      REAL    :: phyt_avg_bkg(:,:,:)
+      REAL    :: tracer_bkg(:,:,:,:)
+      REAL    :: phyto2d_balinc(:,:,:,:)
+      WRITE(*,*) 'asm_phyto3d_bal_medusa: You should not have seen this print! error?'
+   END SUBROUTINE asm_phyto3d_bal_medusa
+#endif
+
+   !!======================================================================
+END MODULE asmphyto3dbal_medusa
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 11989)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90	(revision 11990)
@@ -51,4 +51,6 @@
                                    mld_max, pgrow_avg,                  &
                                    ploss_avg, phyt_avg,                 &
+                                   ploss_avg_3d, phyt_avg_3d,           &
+                                   pgrow_avg_3d,                        &
 # endif
                                    za_sed_c, za_sed_ca, za_sed_fe,      &
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 11989)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90	(revision 11990)
@@ -36,5 +36,6 @@
       USE par_oce,           ONLY: jpi, jpj, jpk
 # if defined key_foam_medusa
-      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
 # else
       USE sms_medusa,        ONLY: jdms
@@ -201,4 +202,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_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 11989)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90	(revision 11990)
@@ -52,4 +52,5 @@
 # if defined key_foam_medusa
                                    pgrow_avg, ploss_avg, phyt_avg,         &
+                                   pgrow_avg_3d, ploss_avg_3d, phyt_avg_3d, &
 # endif
                                    xkphd, xkphn, xkzme, xkzmi,             &
@@ -242,4 +243,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_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 11989)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90	(revision 11990)
@@ -458,4 +458,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
 !!
@@ -537,4 +540,6 @@
 # if defined key_foam_medusa
       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_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 11989)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90	(revision 11990)
@@ -105,4 +105,6 @@
       USE sms_medusa,                 ONLY: hist_pco2, xobs_xco2a,          &
                                             pgrow_avg, ploss_avg,           &
+                                            pgrow_avg_3d, ploss_avg_3d,     &
+                                            phyt_avg_3d,                    &
                                             phyt_avg, mld_max
 # else
Index: /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 11989)
+++ /branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 11990)
@@ -370,4 +370,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
 
@@ -555,4 +566,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
 !!
