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 11958)
+++ 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
