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/asmphyto2dbal_medusa.F90	(revision 11958)
+++ branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto3dbal_medusa.F90	(revision 11990)
@@ -1,3 +1,3 @@
-MODULE asmphyto2dbal_medusa
+MODULE asmphyto3dbal_medusa
    !!======================================================================
    !!                       ***  MODULE asmphyto2dbal_medusa  ***
@@ -34,5 +34,5 @@
    PRIVATE                   
 
-   PUBLIC asm_phyto2d_bal_medusa
+   PUBLIC asm_phyto3d_bal_medusa
 
    ! Default values for biological assimilation parameters
@@ -69,21 +69,10 @@
 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,                    &
+   SUBROUTINE asm_phyto3d_bal_medusa( pinc_chltot,                    &
       &                               pincper,                        &
-      &                               p_maxchlinc, ld_phytobal, pmld, &
+      &                               p_maxchlinc,                    &
       &                               pgrow_avg_bkg, ploss_avg_bkg,   &
-      &                               phyt_avg_bkg, mld_max_bkg,      &
-      &                               tracer_bkg, phyto2d_balinc )
+      &                               phyt_avg_bkg,                   &
+      &                               tracer_bkg, phyto3d_balinc )
       !!---------------------------------------------------------------------------
       !!                    ***  ROUTINE asm_phyto2d_bal_medusa  ***
@@ -101,26 +90,12 @@
       !!---------------------------------------------------------------------------
       !!
-      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(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
-      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)       :: 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) :: phyto2d_balinc ! Balancing increments
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto3d_balinc ! Balancing increments
       !!
       INTEGER                                               :: ji, jj, jk, jn ! Loop counters
@@ -144,306 +119,218 @@
       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(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,jpk,22)    :: diag_fulldepth ! Full-depth 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
-         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
-                  
+                  pinc_chltot(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj,jk), p_maxchlinc ) )
                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
+      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,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
+               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
          
-         ! Propagate through mixed layer
+         ! 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
-               !
-               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
-               !
+
+               ! 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
-
-         ! 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
+      END DO
       
       ! If performing extra tidal mixing in the Indonesian Throughflow,
@@ -453,10 +340,10 @@
          DO jn = 1, jptra
             DO jk = 1, jpk
-               phyto2d_balinc(:,:,jk,jn) = phyto2d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
+               phyto3d_balinc(:,:,jk,jn) = phyto3d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
             END DO
          END DO
       ENDIF
 
-   END SUBROUTINE asm_phyto2d_bal_medusa
+   END SUBROUTINE asm_phyto3d_bal_medusa
 
 #else
@@ -465,47 +352,22 @@
    !!----------------------------------------------------------------------
 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,                    &
+   SUBROUTINE asm_phyto3d_bal_medusa(pinc_chltot,                    &
       &                              pincper,                        &
-      &                              p_maxchlinc, ld_phytobal, pmld, &
+      &                              p_maxchlinc,                    &
       &                              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(:,:)
+      &                              phyt_avg_bkg,                   &
+      &                              tracer_bkg, phyto3d_balinc )
+      REAL    :: pinc_chltot(:,:,:)
       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    :: pgrow_avg_bkg(:,:,:)
+      REAL    :: ploss_avg_bkg(:,:,:)
+      REAL    :: phyt_avg_bkg(:,:,:)
       REAL    :: tracer_bkg(:,:,:,:)
       REAL    :: phyto2d_balinc(:,:,:,:)
-      WRITE(*,*) 'asm_phyto2d_bal_medusa: You should not have seen this print! error?'
-   END SUBROUTINE asm_phyto2d_bal_medusa
+      WRITE(*,*) 'asm_phyto3d_bal_medusa: You should not have seen this print! error?'
+   END SUBROUTINE asm_phyto3d_bal_medusa
 #endif
 
    !!======================================================================
-END MODULE asmphyto2dbal_medusa
+END MODULE asmphyto3dbal_medusa
