Index: branches/UKMO/dev_r5518_GO6_package_FOAMv14_phytobal3d/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FOAMv14_phytobal3d/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90	(revision 12758)
+++ branches/UKMO/dev_r5518_GO6_package_FOAMv14_phytobal3d/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_medusa.F90	(revision 13097)
@@ -1,3 +1,3 @@
-MODULE asmphyto2dbal_medusa
+MODULE asmphytobal_medusa
    !!======================================================================
    !!                       ***  MODULE asmphyto2dbal_medusa  ***
@@ -33,5 +33,5 @@
    PRIVATE                   
 
-   PUBLIC asm_phyto2d_bal_medusa
+   PUBLIC asm_phyto_bal_medusa
 
    ! Default values for biological assimilation parameters
@@ -68,25 +68,26 @@
 CONTAINS
 
-   SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      &
-      &                               pinc_chltot,                    &
-      &                               ld_chldia,                      &
-      &                               pinc_chldia,                    &
-      &                               ld_chlnon,                      &
-      &                               pinc_chlnon,                    &
-      &                               ld_phytot,                      &
-      &                               pinc_phytot,                    &
-      &                               ld_phydia,                      &
-      &                               pinc_phydia,                    &
-      &                               ld_phynon,                      &
-      &                               pinc_phynon,                    &
-      &                               pincper,                        &
-      &                               p_maxchlinc, ld_phytobal, pmld, &
-      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
-      &                               phyt_avg_bkg, mld_max_bkg,      &
-      &                               tracer_bkg, phyto2d_balinc )
+   SUBROUTINE asm_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_phyto2d_bal_medusa  ***
+      !!                    ***  ROUTINE asm_phyto_bal_medusa  ***
       !!
-      !! ** Purpose :   calculate increments to MEDUSA from 2d phytoplankton increments
+      !! ** Purpose :   calculate increments to MEDUSA from phytoplankton increments
       !!
       !! ** Method  :   average up MEDUSA to look like HadOCC
@@ -94,5 +95,5 @@
       !!                separate back out to MEDUSA
       !!
-      !! ** Action  :   populate phyto2d_balinc
+      !! ** Action  :   populate phyto_balinc
       !!
       !! References :   Hemmings et al., 2008, J. Mar. Res.
@@ -100,29 +101,31 @@
       !!---------------------------------------------------------------------------
       !!
-      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,  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
@@ -143,48 +146,52 @@
       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)           :: 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,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
+      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
-         IF ( ld_chltot ) THEN
+         DO jk = 1, kdeps
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
+                  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
-         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
+         END DO
       ENDIF
 
@@ -250,22 +257,24 @@
 
          ! Set background state
-         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin)
-         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd)
-         bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme)
-         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet)
-         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic)
-         bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk)
+         bstate_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(:,:) = 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
+         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
@@ -275,19 +284,62 @@
 
          ! 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 )
+         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
-         phyto2d_balinc(:,:,:,:) = 0.0
+         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
@@ -296,5 +348,5 @@
                   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
+                     & ( pinc_chltot_3d(ji,jj,jkinc) /= 0.0 ) ) THEN
                      IF ( ld_chltot ) THEN
                         ! Phytoplankton nitrogen split up based on existing ratios
@@ -305,6 +357,6 @@
                      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)
+                        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
 
@@ -318,9 +370,9 @@
                      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
+                     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
 
@@ -331,23 +383,23 @@
                      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
+                     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
-                  phyto2d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1))
+                  phyto_balinc(ji,jj,jk,jpdin) = outincs_3d(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))
+                  phyto_balinc(ji,jj,jk,jpdet) = outincs_3d(ji,jj,jk,i_tracer(4))
 
                   ! DIC straight from balancing scheme
-                  phyto2d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5))
+                  phyto_balinc(ji,jj,jk,jpdic) = outincs_3d(ji,jj,jk,i_tracer(5))
 
                   ! Alkalinity straight from balancing scheme
-                  phyto2d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6))
+                  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) - phyto2d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
-                     phyto2d_balinc(ji,jj,jk,jpsil) = phyto2d_balinc(ji,jj,jk,jppds) * (-1.0)
+                  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
 
@@ -355,10 +407,10 @@
                   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
+                     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
-                  phyto2d_balinc(ji,jj,jk,jpfer) = 0.0
-                  phyto2d_balinc(ji,jj,jk,jpoxy) = 0.0
+                  phyto_balinc(ji,jj,jk,jpfer) = 0.0
+                  phyto_balinc(ji,jj,jk,jpoxy) = 0.0
                   
                END DO
@@ -369,78 +421,82 @@
       
          ! Initialise individual chlorophyll increments to zero
-         phyto2d_balinc(:,:,:,jpchn) = 0.0
-         phyto2d_balinc(:,:,:,jpchd) = 0.0
+         phyto_balinc(:,:,:,jpchn) = 0.0
+         phyto_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
+         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
          
-         ! 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
+         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
-               !
-               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
+            END DO
+         ENDIF
 
          ! 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
+         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
@@ -452,10 +508,10 @@
          DO jn = 1, jptra
             DO jk = 1, jpk
-               phyto2d_balinc(:,:,jk,jn) = phyto2d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
+               phyto_balinc(:,:,jk,jn) = phyto_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
             END DO
          END DO
       ENDIF
 
-   END SUBROUTINE asm_phyto2d_bal_medusa
+   END SUBROUTINE asm_phyto_bal_medusa
 
 #else
@@ -464,47 +520,49 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      &
-      &                              pinc_chltot,                    &
-      &                              ld_chldia,                      &
-      &                              pinc_chldia,                    &
-      &                              ld_chlnon,                      &
-      &                              pinc_chlnon,                    &
-      &                              ld_phytot,                      &
-      &                              pinc_phytot,                    &
-      &                              ld_phydia,                      &
-      &                              pinc_phydia,                    &
-      &                              ld_phynon,                      &
-      &                              pinc_phynon,                    &
-      &                              pincper,                        &
-      &                              p_maxchlinc, ld_phytobal, pmld, &
-      &                              pgrow_avg_bkg, ploss_avg_bkg,   &
-      &                              phyt_avg_bkg, mld_max_bkg,      &
-      &                              tracer_bkg, phyto2d_balinc )
+   SUBROUTINE asm_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(:,:)
+      REAL    :: pinc_chltot_3d(:,:,:)
       LOGICAL :: ld_chldia
-      REAL    :: pinc_chldia(:,:)
+      REAL    :: pinc_chldia_3d(:,:,:)
       LOGICAL :: ld_chlnon
-      REAL    :: pinc_chlnon(:,:)
+      REAL    :: pinc_chlnon_3d(:,:,:)
       LOGICAL :: ld_phytot
-      REAL    :: pinc_phytot(:,:)
+      REAL    :: pinc_phytot_3d(:,:,:)
       LOGICAL :: ld_phydia
-      REAL    :: pinc_phydia(:,:)
+      REAL    :: pinc_phydia_3d(:,:,:)
       LOGICAL :: ld_phynon
-      REAL    :: pinc_phynon(:,:)
+      REAL    :: pinc_phynon_3d(:,:,:)
       REAL    :: pincper
       REAL    :: p_maxchlinc
       LOGICAL :: ld_phytobal
       REAL    :: pmld(:,:)
-      REAL    :: pgrow_avg_bkg(:,:)
-      REAL    :: ploss_avg_bkg(:,:)
-      REAL    :: phyt_avg_bkg(:,:)
+      REAL    :: pgrow_avg_bkg_3d(:,:,:)
+      REAL    :: ploss_avg_bkg_3d(:,:,:)
+      REAL    :: phyt_avg_bkg_3d(:,:,:)
       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
+      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 asmphyto2dbal_medusa
+END MODULE asmphytobal_medusa
