diff --git a/src/accessom_coupler/mom_oasis3_interface.F90 b/src/accessom_coupler/mom_oasis3_interface.F90
index 77c729b4e2..939ecf8319 100644
--- a/src/accessom_coupler/mom_oasis3_interface.F90
+++ b/src/accessom_coupler/mom_oasis3_interface.F90
@@ -155,9 +155,9 @@ module mom_oasis3_interface_mod
 integer iisc,iiec,jjsc,jjec
 integer iisd,iied,jjsd,jjed
 
-integer, parameter :: max_fields_in=23
+integer, parameter :: max_fields_in=25
 
-integer, parameter :: max_fields_out=8
+integer, parameter :: max_fields_out=10
 
 integer, dimension(max_fields_in)  :: id_var_in  ! ID for fields to be rcvd
 integer, dimension(max_fields_out) :: id_var_out ! ID for fields to be sent
@@ -375,6 +375,10 @@ subroutine coupler_init(Dom, Time, Time_step_coupled, dt_cpld, Run_len, &
   mom_name_read(21)='licefw'  ! Water flux from land ice
   mom_name_read(22)='liceht'  ! Heat flux from land ice
   mom_name_read(23)='wnd_io'  !
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  mom_name_read(24)='iof_nit'  !
+  mom_name_read(25)='iof_alg'  !
+#endif
 
   !ocn ==> ice
   mom_name_write(:)=''
@@ -387,7 +391,10 @@ subroutine coupler_init(Dom, Time, Time_step_coupled, dt_cpld, Run_len, &
   mom_name_write(6)='frazil'
   mom_name_write(7)='dssldx'
   mom_name_write(8)='dssldy'
-
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  mom_name_write(9)='n_surf'
+  mom_name_write(10)='alg_surf'
+#endif
 
   fmatch = .false.
   do jf = 1,num_fields_in
@@ -597,6 +604,12 @@ subroutine into_coupler(step, Ocean_sfc, Time, before_ocean_update)
     vtmp(iisd:iied,jjsd:jjed) = Ocean_sfc%gradient(iisd:iied,jjsd:jjed,1)
   case('dssldy') 
     vtmp(iisd:iied,jjsd:jjed) = Ocean_sfc%gradient(iisd:iied,jjsd:jjed,2)
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  case('n_surf')
+    vtmp(iisd:iied,jjsd:jjed) = Ocean_sfc%n_surf(iisd:iied,jjsd:jjed)
+  case('alg_surf')
+    vtmp(iisd:iied,jjsd:jjed) = Ocean_sfc%alg_surf(iisd:iied,jjsd:jjed)
+#endif
   case DEFAULT
   call mpp_error(FATAL,&
       '==>Error from into_coupler: Unknown quantity.')
@@ -746,6 +759,12 @@ subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time)
      Ice_ocean_boundary%fprec(iisc:iiec,jjsc:jjec) =  vwork(iisc:iiec,jjsc:jjec)
   case('aice')
      Ice_ocean_boundary%aice(iisc:iiec,jjsc:jjec) =  vwork(iisc:iiec,jjsc:jjec)
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  case('iof_nit')
+     Ice_ocean_boundary%iof_nit(iisc:iiec,jjsc:jjec) =  vwork(iisc:iiec,jjsc:jjec)
+  case('iof_alg')
+     Ice_ocean_boundary%iof_alg(iisc:iiec,jjsc:jjec) =  vwork(iisc:iiec,jjsc:jjec)
+#endif
   case('mh_flux')
      Ice_ocean_boundary%mh_flux(iisc:iiec,jjsc:jjec) =  vwork(iisc:iiec,jjsc:jjec)
   case('wfimelt')
@@ -820,6 +839,10 @@ subroutine write_coupler_restart(step,Ocean_sfc,write_restart)
           case('dssldx'); vtmp = Ocean_sfc%gradient(iisd:iied,jjsd:jjed,1); fld_ice='sslx_i'
           case('dssldy'); vtmp = Ocean_sfc%gradient(iisd:iied,jjsd:jjed,2); fld_ice='ssly_i'
           case('frazil'); vtmp = Ocean_sfc%frazil(iisd:iied,jjsd:jjed); fld_ice='pfmice_i'
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+          case('n_surf'); vtmp = Ocean_sfc%n_surf(iisd:iied,jjsd:jjed); fld_ice='ssn_i'
+          case('alg_surf'); vtmp = Ocean_sfc%alg_surf(iisd:iied,jjsd:jjed); fld_ice='ssalg_i'
+#endif
         end select
 
         if (parallel_coupling) then
diff --git a/src/accessom_coupler/ocean_solo.F90 b/src/accessom_coupler/ocean_solo.F90
index 52d2f474f4..838759cde0 100644
--- a/src/accessom_coupler/ocean_solo.F90
+++ b/src/accessom_coupler/ocean_solo.F90
@@ -413,7 +413,10 @@ program main
              Ice_ocean_boundary% licefw(isc:iec,jsc:jec),           &
              Ice_ocean_boundary% liceht(isc:iec,jsc:jec),           &
              Ice_ocean_boundary%wnd(isc:iec,jsc:jec))
-
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  allocate ( Ice_ocean_boundary%iof_nit(isc:iec,jsc:jec),           &
+             Ice_ocean_boundary%iof_alg(isc:iec,jsc:jec))
+#endif
   Ice_ocean_boundary%u_flux          = 0.0
   Ice_ocean_boundary%v_flux          = 0.0
   Ice_ocean_boundary%t_flux          = 0.0
@@ -436,7 +439,10 @@ program main
   Ice_ocean_boundary%licefw          = 0.0
   Ice_ocean_boundary%liceht          = 0.0
   Ice_ocean_boundary%wnd             = 0.0
-
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  Ice_ocean_boundary%iof_nit         = 0.0
+  Ice_ocean_boundary%iof_alg         = 0.0
+#endif
   coupler_init_clock = mpp_clock_id('OASIS init', grain=CLOCK_COMPONENT)
   call mpp_clock_begin(coupler_init_clock)
   call external_coupler_sbc_init(Ocean_sfc%domain, dt_cpld, Run_len, &
@@ -702,7 +708,10 @@ subroutine write_boundary_chksums(Ice_ocean_boundary)
     call write_chksum_2d('Ice_ocean_boundary%wfimelt', Ice_ocean_boundary%wfimelt)
     call write_chksum_2d('Ice_ocean_boundary%wfiform', Ice_ocean_boundary%wfiform)
     call write_chksum_2d('Ice_ocean_boundary%wnd', Ice_ocean_boundary%wnd)
-
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+    call write_chksum_2d('Ice_ocean_boundary%iof_nit', Ice_ocean_boundary%iof_nit)
+    call write_chksum_2d('Ice_ocean_boundary%iof_alg', Ice_ocean_boundary%iof_alg)
+#endif
 end subroutine write_boundary_chksums
 
 end program main
diff --git a/src/mom5/ocean_core/ocean_sbc.F90 b/src/mom5/ocean_core/ocean_sbc.F90
index fcd790d19a..d0ac92f0a5 100644
--- a/src/mom5/ocean_core/ocean_sbc.F90
+++ b/src/mom5/ocean_core/ocean_sbc.F90
@@ -535,7 +535,7 @@ module ocean_sbc_mod
 use ocean_tracer_util_mod,    only: diagnose_3d_rho
 
 #if defined(CSIRO_BGC)
-use csiro_bgc_mod,            only: csiro_bgc_virtual_fluxes, do_csiro_bgc
+use csiro_bgc_mod,            only: csiro_bgc_virtual_fluxes, do_csiro_bgc,ind_no3,ind_phy
 #endif
 implicit none
 
@@ -681,7 +681,10 @@ module ocean_sbc_mod
 integer :: id_mh_flux            =-1
 integer :: id_atm_co2            =-1
 #endif
-
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+integer :: id_iof_nit               =-1
+integer :: id_iof_alg               =-1
+#endif
 
 ! ids for sea level forcing fields 
 real    :: cellarea_r
@@ -821,7 +824,7 @@ module ocean_sbc_mod
 #endif
 #if defined(ACCESS_CM) || defined(ACCESS_OM)
 real, allocatable, dimension(:,:,:) :: sslope
-real, allocatable, dimension(:,:) :: aice
+real, allocatable, dimension(:,:) :: aice, iof_nit, iof_alg
 #endif
 #if defined(ACCESS_CM)
 real, allocatable, dimension(:,:) :: co2flux
@@ -1116,12 +1119,18 @@ subroutine ocean_sbc_init(Grid, Domain, Time, T_prog, T_diag, &
   allocate ( Ocean_sfc%gradient (isc_bnd:iec_bnd,jsc_bnd:jec_bnd,2))
   allocate ( sslope(isc:iec, jsc:jec, 2) )
   allocate ( aice(isd:ied, jsd:jed) )
+#endif
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  allocate ( iof_nit(isd:ied, jsd:jed) )
+  allocate ( iof_alg(isd:ied, jsd:jed) )
+  allocate ( Ocean_sfc%n_surf (isc_bnd:iec_bnd,jsc_bnd:jec_bnd))
+  allocate ( Ocean_sfc%alg_surf (isc_bnd:iec_bnd,jsc_bnd:jec_bnd))
+#endif
 #if defined(ACCESS_CM)
   allocate ( Ocean_sfc%co2    (isc_bnd:iec_bnd,jsc_bnd:jec_bnd), &
              Ocean_sfc%co2flux (isc_bnd:iec_bnd,jsc_bnd:jec_bnd)) 
   allocate ( co2flux(isd:ied,jsd:jed),ocn_co2(isd:ied,jsd:jed))
   allocate ( atm_co2(isd:ied,jsd:jed))
-#endif
 #endif
 
   Ocean_sfc%t_surf  = 0.0  ! time averaged sst (Kelvin) passed to atmosphere/ice model
@@ -1136,6 +1145,12 @@ subroutine ocean_sbc_init(Grid, Domain, Time, T_prog, T_diag, &
   sslope = 0.0
   aice = 0.0
 #endif
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  iof_nit = 0.0
+  iof_alg = 0.0
+  Ocean_sfc%n_surf  = 0.0 
+  Ocean_sfc%alg_surf  = 0.0 
+#endif
 #if defined(ACCESS_CM)
   Ocean_sfc%co2       = 0.0 
   Ocean_sfc%co2flux   = 0.0 
@@ -1992,6 +2007,16 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog)
        Time%model_time, 'fraction of surface area covered with ice', 'm^2/m^2' ,  &
        missing_value=missing_value,range=(/-1.e1,1.e1/),                      &
        standard_name='areal_ice_concentration' )
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  id_iof_nit = register_diag_field('ocean_model','iof_nit', Grd%tracer_axes(1:2),&
+       Time%model_time, 'ice-ocean flux of nitrate', 'mmol/m^2/s^1' ,  &
+       missing_value=missing_value,range=(/-1.e1,1.e1/),                      &
+       standard_name='ice_ocean_nitrate_flux' )
+  id_iof_alg = register_diag_field('ocean_model','iof_alg', Grd%tracer_axes(1:2),&
+       Time%model_time, 'ice-ocean flux of algae', 'mmol/m^2/s^1' ,  &
+       missing_value=missing_value,range=(/-1.e1,1.e1/),                      &
+       standard_name='ice_ocean_algal_flux' )
+#endif
   id_wnd = register_diag_field('ocean_model','wnd', Grd%tracer_axes(1:2),&
        Time%model_time, 'Wind speed', 'm/s' ,  &
        missing_value=missing_value,range=(/-1.e3,1.e3/),                      &
@@ -2819,6 +2844,15 @@ subroutine initialize_ocean_sfc(Time, Thickness, T_prog, T_diag, Velocity, Ocean
 #endif
   end where
 
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  if (ind_no3 > 0) then
+   where (Grd%tmask(isc:iec,jsc:jec,1) == 1.0)
+      Ocean_sfc%n_surf(isc_bnd:iec_bnd,jsc_bnd:jec_bnd)  = T_prog(ind_no3)%field(isc:iec,jsc:jec,1,taup1)
+      Ocean_sfc%alg_surf(isc_bnd:iec_bnd,jsc_bnd:jec_bnd)  = T_prog(ind_phy)%field(isc:iec,jsc:jec,1,taup1)
+   end where
+  end if
+#endif
+
   ! when enabled, use FAFMIP redistributed heat tracer for sst
   if(index_redist_heat > 0) then
     write(stdoutunit,*) &
@@ -2851,6 +2885,10 @@ subroutine initialize_ocean_sfc(Time, Thickness, T_prog, T_diag, Velocity, Ocean
   id_field = register_restart_field(Sfc_restart, filename, 'v_surf', Ocean_sfc%v_surf,Ocean_sfc%Domain)
   id_field = register_restart_field(Sfc_restart, filename, 'sea_lev',Ocean_sfc%sea_lev,Ocean_sfc%Domain)
   id_field = register_restart_field(Sfc_restart, filename, 'frazil', Ocean_sfc%frazil,Ocean_sfc%Domain)
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  id_field = register_restart_field(Sfc_restart, filename, 'n_surf', Ocean_sfc%n_surf,Ocean_sfc%Domain)
+  id_field = register_restart_field(Sfc_restart, filename, 'alg_surf', Ocean_sfc%alg_surf,Ocean_sfc%Domain)
+#endif
 #if defined(ACCESS_CM)
 !RASF Make these optional so we don't break existing runs.
   id_field = register_restart_field(Sfc_restart, filename, 'co2flux',Ocean_sfc%co2flux,Ocean_sfc%Domain, mandatory=.false.)
@@ -2941,6 +2979,12 @@ subroutine sum_ocean_sfc(Time, Thickness, T_prog, T_diag, Dens, Velocity, Ocean_
 #if defined(ACCESS_CM) || defined(ACCESS_OM)
             Ocean_sfc%gradient(i,j,:) = Ocean_sfc%gradient(i,j,:) + sslope(ii,jj,:)  
 #endif
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+            if (ind_no3 > 0) then
+             Ocean_sfc%n_surf(i,j)  = Ocean_sfc%n_surf(i,j)  + T_prog(ind_no3)%field(ii,jj,1,taup1)
+             Ocean_sfc%alg_surf(i,j)  = Ocean_sfc%alg_surf(i,j)  + T_prog(ind_phy)%field(ii,jj,1,taup1)
+            end if
+#endif
 #if defined(ACCESS_CM)
             Ocean_sfc%co2flux(i,j) = Ocean_sfc%co2flux(i,j) + co2flux(ii,jj)
             Ocean_sfc%co2(i,j)     = Ocean_sfc%co2(i,j) + ocn_co2(ii,jj)
@@ -2981,6 +3025,12 @@ subroutine sum_ocean_sfc(Time, Thickness, T_prog, T_diag, Dens, Velocity, Ocean_
             Ocean_sfc%v_surf(i,j)  = Ocean_sfc%v_surf(i,j)  + Grd%umask(ii,jj,1)*onehalf &
                                      *(Velocity%u(ii,jj,1,2,taup1) + Velocity%u(ii,jj+1,1,2,taup1))
             Ocean_sfc%sea_lev(i,j) = Ocean_sfc%sea_lev(i,j) + Thickness%sea_lev(ii,jj)
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+            if (ind_no3 > 0) then
+             Ocean_sfc%n_surf(i,j)  = Ocean_sfc%n_surf(i,j)  + T_prog(ind_no3)%field(ii,jj,1,taup1)
+             Ocean_sfc%alg_surf(i,j)  = Ocean_sfc%alg_surf(i,j)  + T_prog(ind_phy)%field(ii,jj,1,taup1)
+            end if
+#endif         
          enddo
       enddo
 
@@ -3044,6 +3094,10 @@ subroutine zero_ocean_sfc(Ocean_sfc)
 #if defined(ACCESS_CM) || defined(ACCESS_OM)
         Ocean_sfc%gradient(i,j,:)= 0.0
 #endif
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+        Ocean_sfc%n_surf(i,j)= 0.0
+        Ocean_sfc%alg_surf(i,j)= 0.0
+#endif
 #if defined(ACCESS_CM)
         Ocean_sfc%co2flux(i,j) = 0.0
         Ocean_sfc%co2(i,j)     = 0.0
@@ -3125,6 +3179,10 @@ subroutine avg_ocean_sfc(Time, Thickness, T_prog, T_diag, Velocity, Ocean_sfc)
 #if defined(ACCESS_CM) || defined(ACCESS_OM)
            Ocean_sfc%gradient(i,j,:) = Ocean_sfc%gradient(i,j,:)*divid
 #endif
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+           Ocean_sfc%n_surf(i,j)  = Ocean_sfc%n_surf(i,j)*divid
+           Ocean_sfc%alg_surf(i,j)  = Ocean_sfc%alg_surf(i,j)*divid
+#endif
 #if defined(ACCESS_CM)
            Ocean_sfc%co2flux(i,j) = Ocean_sfc%co2flux(i,j)*divid
            Ocean_sfc%co2(i,j)     = Ocean_sfc%co2(i,j)*divid
@@ -4351,6 +4409,10 @@ subroutine get_ocean_sbc(Time, Ice_ocean_boundary, Thickness, Dens, Ext_mode, T_
         ii = i + i_shift
         jj = j + j_shift
         aice(ii,jj) = Ice_ocean_boundary%aice(i,j)*Grd%tmask(ii,jj,1)
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+        iof_nit(ii,jj) = Ice_ocean_boundary%iof_nit(i,j)*Grd%tmask(ii,jj,1)
+        iof_alg(ii,jj) = Ice_ocean_boundary%iof_alg(i,j)*Grd%tmask(ii,jj,1)
+#endif
      enddo
   enddo
 #endif
@@ -4365,7 +4427,7 @@ subroutine get_ocean_sbc(Time, Ice_ocean_boundary, Thickness, Dens, Ext_mode, T_
 #elif defined(ACCESS_OM) && defined(CSIRO_BGC)
 ! Do not pass co2flux, ocn_co2 or atm_co2
   call ocean_tpm_sbc(Dom, Grd, T_prog(:), Time, Ice_ocean_boundary%fluxes, runoff, &
-                     isc_bnd, iec_bnd, jsc_bnd, jec_bnd,aice=aice, wnd=Velocity%u10, &
+                     isc_bnd, iec_bnd, jsc_bnd, jec_bnd,aice=aice, iof_nit=iof_nit, iof_alg=iof_alg, wnd=Velocity%u10, &
                      use_waterflux=use_waterflux, salt_restore_as_salt_flux=salt_restore_as_salt_flux)
 #else 
   call ocean_tpm_sbc(Dom, Grd, T_prog(:), Time, Ice_ocean_boundary%fluxes, runoff, &
@@ -6022,6 +6084,14 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou
   if (id_aice > 0) used = send_data(id_aice, aice(:,:),    &
                  Time%model_time, rmask=Grd%tmask(:,:,1),  &
                  is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+  if (id_iof_nit > 0) used = send_data(id_iof_nit, iof_nit(:,:),    &
+                 Time%model_time, rmask=Grd%tmask(:,:,1),  &
+                 is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
+  if (id_iof_alg > 0) used = send_data(id_iof_alg, iof_alg(:,:),    &
+                 Time%model_time, rmask=Grd%tmask(:,:,1),  &
+                 is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
+#endif
   if (id_wnd > 0) used = send_data(id_wnd, Velocity%u10(:,:),    &
                  Time%model_time, rmask=Grd%tmask(:,:,1),  &
                  is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
diff --git a/src/mom5/ocean_core/ocean_types.F90 b/src/mom5/ocean_core/ocean_types.F90
index 44acce0d68..2ed795e1b6 100644
--- a/src/mom5/ocean_core/ocean_types.F90
+++ b/src/mom5/ocean_core/ocean_types.F90
@@ -1244,6 +1244,10 @@ module ocean_types_mod
      real, pointer, dimension(:,:) :: co2              =>NULL() ! co2
 #endif
      real, pointer, dimension(:,:) :: wnd              =>NULL() ! wind speed
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+     real, pointer, dimension(:,:) :: iof_nit              =>NULL() ! ice-ocean flux of nitrate
+     real, pointer, dimension(:,:) :: iof_alg              =>NULL() ! ice-ocean flux of algae
+#endif
      integer :: xtype                                          ! REGRID, REDIST or DIRECT
 
      type(coupler_2d_bc_type)      :: fluxes                   ! array of fields used for additional tracers
@@ -1266,6 +1270,10 @@ module ocean_types_mod
 #if defined(ACCESS_CM) || defined(ACCESS_OM)
      real, pointer, dimension(:,:,:)  :: gradient =>NULL() ! x/y slopes of sea surface.
 #endif
+#if defined(ACCESS_OM) && defined(CSIRO_BGC)
+     real, pointer, dimension(:,:)    :: n_surf =>NULL() ! sea surface nitrate (mmol m-3)
+     real, pointer, dimension(:,:)    :: alg_surf =>NULL() ! sea surface algae (mmol m-3)
+#endif
 #if defined(ACCESS_CM)
      real, pointer, dimension(:,:)    :: co2     =>NULL() ! co2 ( )
      real, pointer, dimension(:,:)    :: co2flux =>NULL() ! co2 flux ()
diff --git a/src/mom5/ocean_csiro_bgc/bio_v3.inc b/src/mom5/ocean_csiro_bgc/bio_v3.inc
index 59558c5a3c..ef65b96234 100755
--- a/src/mom5/ocean_csiro_bgc/bio_v3.inc
+++ b/src/mom5/ocean_csiro_bgc/bio_v3.inc
@@ -141,6 +141,29 @@ do n = 1, instances  !{
   pprod_gross(:,:,:) = 0.0
   zprod_gross(:,:,:) = 0.0
   light_limit(:,:) = 0.0
+  radbio3d(:,:,:) = 0.0
+  npp3d(:,:,:) = 0.0
+  adic_intmld(:,:) = 0.0
+  dic_intmld(:,:) = 0.0
+  o2_intmld(:,:) = 0.0
+  no3_intmld(:,:) = 0.0
+  fe_intmld(:,:) = 0.0
+  phy_intmld(:,:) = 0.0
+  det_intmld(:,:) = 0.0
+  pprod_gross_intmld(:,:) = 0.0
+  npp_intmld(:,:) = 0.0
+  radbio_intmld(:,:) = 0.0
+  adic_int100(:,:) = 0.0
+  dic_int100(:,:) = 0.0
+  o2_int100(:,:) = 0.0
+  no3_int100(:,:) = 0.0
+  fe_int100(:,:) = 0.0
+  phy_int100(:,:) = 0.0
+  det_int100(:,:) = 0.0
+  pprod_gross_int100(:,:) = 0.0
+  npp_int100(:,:) = 0.0
+  radbio_int100(:,:) = 0.0
+
 
 ! possible tracers in model
   ind_dic = biotic(n)%ind_bgc(id_dic)
@@ -223,6 +246,8 @@ do n = 1, instances  !{
 !          radbio = max(0.0,parbio(i,j)*sw_zt)
 ! well, now doing an experiment without ice, so assumption no longer true.  mac, oct12. 
 
+          radbio3d(i,j,k) = radbio
+
           vpbio(i,k) = abio(i,j)*bbio(i,j)**                                  &
                   (cbio(i,j)*t_prog(index_temp)%field(i,j,k,time%tau))
 
@@ -319,6 +344,8 @@ do n = 1, instances  !{
 
 ! Estimate primary productivity from phytoplankton growth
           pprod_gross(i,j,k) = pprod_gross(i,j,k) + dtsb*f11
+! Net primary productivity (gross PP minus linear mortality)
+	  npp3d(i,j,k) = npp3d(i,j,k) + dtsb*(f11 - f21)
 
 !chd   -- zooplankton equation
           biotr(i,k,id_zoo)  = biotr(i,k,id_zoo)  + dtsb * (                     &
@@ -393,7 +420,33 @@ do n = 1, instances  !{
 
         pprod_gross(i,j,k)=rdtts*pprod_gross(i,j,k)*grid%tmask(i,j,k)
         zprod_gross(i,j,k)=rdtts*zprod_gross(i,j,k)*grid%tmask(i,j,k)
+        npp3d(i,j,k)=rdtts*npp3d(i,j,k)*grid%tmask(i,j,k)
+
+          if (Grid%zw(k) .le. mld(i,j)) then
+            adic_intmld(i,j)= adic_intmld(i,j)+ t_prog(ind_adic)%field(i,j,k,time%taum1)* thickness%dzt(i,j,k)
+            dic_intmld(i,j) = dic_intmld(i,j) + t_prog(ind_dic)%field(i,j,k,time%taum1) * thickness%dzt(i,j,k)
+            o2_intmld(i,j)  = o2_intmld(i,j)  + t_prog(ind_o2)%field(i,j,k,time%taum1)  * thickness%dzt(i,j,k)
+            no3_intmld(i,j) = no3_intmld(i,j) + t_prog(ind_no3)%field(i,j,k,time%taum1) * thickness%dzt(i,j,k)
+            fe_intmld(i,j)  = fe_intmld(i,j)  + t_prog(ind_fe)%field(i,j,k,time%taum1)  * thickness%dzt(i,j,k)
+            phy_intmld(i,j) = phy_intmld(i,j) + t_prog(ind_phy)%field(i,j,k,time%taum1) * thickness%dzt(i,j,k)
+            det_intmld(i,j) = det_intmld(i,j) + t_prog(ind_det)%field(i,j,k,time%taum1) * thickness%dzt(i,j,k)
+            pprod_gross_intmld(i,j) = pprod_gross_intmld(i,j) + pprod_gross(i,j,k) * thickness%dzt(i,j,k)
+            npp_intmld(i,j) = npp_intmld(i,j) + npp3d(i,j,k) * thickness%dzt(i,j,k)
+            radbio_intmld(i,j) = radbio_intmld(i,j) + radbio3d(i,j,k) * thickness%dzt(i,j,k)
+          endif
 
+          if (Grid%zw(k) .le. 100) then
+            adic_int100(i,j)= adic_int100(i,j)+ t_prog(ind_adic)%field(i,j,k,time%taum1)* thickness%dzt(i,j,k)
+            dic_int100(i,j) = dic_int100(i,j) + t_prog(ind_dic)%field(i,j,k,time%taum1) * thickness%dzt(i,j,k)
+            o2_int100(i,j)  = o2_int100(i,j)  + t_prog(ind_o2)%field(i,j,k,time%taum1)  * thickness%dzt(i,j,k)
+            no3_int100(i,j) = no3_int100(i,j) + t_prog(ind_no3)%field(i,j,k,time%taum1) * thickness%dzt(i,j,k)
+            fe_int100(i,j)  = fe_int100(i,j)  + t_prog(ind_fe)%field(i,j,k,time%taum1)  * thickness%dzt(i,j,k)
+            phy_int100(i,j) = phy_int100(i,j) + t_prog(ind_phy)%field(i,j,k,time%taum1) * thickness%dzt(i,j,k)
+            det_int100(i,j) = det_int100(i,j) + t_prog(ind_det)%field(i,j,k,time%taum1) * thickness%dzt(i,j,k)
+            pprod_gross_int100(i,j) = pprod_gross_int100(i,j) + pprod_gross(i,j,k) * thickness%dzt(i,j,k)
+            npp_int100(i,j) = npp_int100(i,j) + npp3d(i,j,k) * thickness%dzt(i,j,k)
+            radbio_int100(i,j) = radbio_int100(i,j) + radbio3d(i,j,k) * thickness%dzt(i,j,k)
+          endif
 
          enddo  !} i
         enddo  !} k
@@ -407,6 +460,7 @@ do n = 1, instances  !{
       enddo
      enddo
 
+
 !RASF no flux boundary conditions
      do i=isc,iec
       adv_fb(i,1)               = 0.0
diff --git a/src/mom5/ocean_csiro_bgc/csiro_bgc.F90 b/src/mom5/ocean_csiro_bgc/csiro_bgc.F90
index af4e7178fa..ef88d8f5d6 100755
--- a/src/mom5/ocean_csiro_bgc/csiro_bgc.F90
+++ b/src/mom5/ocean_csiro_bgc/csiro_bgc.F90
@@ -221,8 +221,17 @@ module  csiro_bgc_mod  !{
 integer :: id_po4, id_dic, id_alk, id_o2, id_no3, id_phy, id_det, id_zoo &
       , id_caco3, id_adic, id_fe, id_caco3_sediment, id_det_sediment
 ! internal pointer to make reading the code easier
-integer :: ind_po4, ind_dic, ind_alk, ind_o2, ind_no3, ind_phy, ind_det, ind_zoo &
-      , ind_caco3, ind_adic, ind_fe
+integer,public :: ind_po4 = -1
+integer,public :: ind_dic = -1 
+integer,public :: ind_alk = -1
+integer,public :: ind_o2 = -1
+integer,public :: ind_no3 = -1
+integer,public :: ind_phy = -1
+integer,public :: ind_det = -1
+integer,public :: ind_zoo = -1
+integer,public :: ind_caco3 = -1
+integer,public :: ind_adic = -1
+integer,public :: ind_fe = -1
 character*6  :: qbio_model
 integer      :: bio_version    ! version of the bgc module to use
 logical      :: zero_floor     ! apply hard floor to bgc tracers 
@@ -238,6 +247,32 @@ module  csiro_bgc_mod  !{
 character*32                            :: atmpress_name    
 real, allocatable, dimension(:,:)       :: fice_t
 integer                                 :: id_light_limit = -1
+integer                                 :: id_adic_intmld = -1
+integer                                 :: id_dic_intmld = -1
+integer                                 :: id_o2_intmld = -1
+integer                                 :: id_no3_intmld = -1
+integer                                 :: id_fe_intmld = -1
+integer                                 :: id_phy_intmld = -1
+integer                                 :: id_det_intmld = -1
+integer                                 :: id_pprod_gross_intmld = -1
+integer                                 :: id_npp_intmld = -1
+integer                                 :: id_radbio_intmld = -1
+integer                                 :: id_adic_int100 = -1
+integer                                 :: id_dic_int100 = -1
+integer                                 :: id_o2_int100 = -1
+integer                                 :: id_no3_int100 = -1
+integer                                 :: id_fe_int100 = -1
+integer                                 :: id_phy_int100 = -1
+integer                                 :: id_det_int100 = -1
+integer                                 :: id_pprod_gross_int100 = -1
+integer                                 :: id_npp_int100 = -1
+integer                                 :: id_radbio_int100 = -1
+integer                                 :: id_radbio1 = -1
+integer                                 :: id_radbio3d = -1
+integer                                 :: id_wdet100 = -1
+integer                                 :: id_npp1 = -1
+integer                                 :: id_npp2d = -1
+integer                                 :: id_npp3d = -1
 integer                                 :: id_pprod_gross = -1
 integer                                 :: id_pprod_gross_2d = -1
 integer                                 :: id_zprod_gross = -1
@@ -303,6 +338,14 @@ module  csiro_bgc_mod  !{
 real, allocatable, dimension(:) :: fmin_pic
 real, allocatable, dimension(:,:,:) :: biotr
 real, allocatable, dimension(:,:) :: light_limit
+real, allocatable, dimension(:,:) :: adic_intmld,dic_intmld,o2_intmld,no3_intmld,fe_intmld,phy_intmld,det_intmld
+real, allocatable, dimension(:,:) :: adic_int100,dic_int100,o2_int100,no3_int100,fe_int100,phy_int100,det_int100
+real, allocatable, dimension(:,:) :: pprod_gross_intmld,npp_intmld,radbio_intmld
+real, allocatable, dimension(:,:) :: pprod_gross_int100,npp_int100,radbio_int100
+real, allocatable, dimension(:,:,:) :: radbio3d
+real, allocatable, dimension(:,:) :: wdet100
+real, allocatable, dimension(:,:) :: npp2d
+real, allocatable, dimension(:,:,:) :: npp3d
 real, allocatable, dimension(:,:,:) :: pprod_gross
 real, allocatable, dimension(:,:) :: pprod_gross_2d
 real, allocatable, dimension(:,:,:) :: zprod_gross
@@ -451,6 +494,30 @@ subroutine allocate_arrays (isc, iec, jsc, jec, isd, ied, jsd, jed, nk)  !{
 allocate( ray(nk) )
 allocate( biotr(isc:iec,nk,ntr_bgc) )
 allocate( light_limit(isc:iec,jsc:jec) )
+allocate( adic_intmld(isc:iec,jsc:jec) )
+allocate( dic_intmld(isc:iec,jsc:jec) )
+allocate( o2_intmld(isc:iec,jsc:jec) )
+allocate( fe_intmld(isc:iec,jsc:jec) )
+allocate( no3_intmld(isc:iec,jsc:jec) )
+allocate( phy_intmld(isc:iec,jsc:jec) )
+allocate( det_intmld(isc:iec,jsc:jec) )
+allocate( pprod_gross_intmld(isc:iec,jsc:jec) )
+allocate( npp_intmld(isc:iec,jsc:jec) )
+allocate( radbio_intmld(isc:iec,jsc:jec) )
+allocate( adic_int100(isc:iec,jsc:jec) )
+allocate( dic_int100(isc:iec,jsc:jec) )
+allocate( o2_int100(isc:iec,jsc:jec) )
+allocate( fe_int100(isc:iec,jsc:jec) )
+allocate( no3_int100(isc:iec,jsc:jec) )
+allocate( phy_int100(isc:iec,jsc:jec) )
+allocate( det_int100(isc:iec,jsc:jec) )
+allocate( pprod_gross_int100(isc:iec,jsc:jec) )
+allocate( npp_int100(isc:iec,jsc:jec) )
+allocate( radbio_int100(isc:iec,jsc:jec) )
+allocate( radbio3d(isc:iec,jsc:jec,nk) )
+allocate( wdet100(isc:iec,jsc:jec) )
+allocate( npp2d(isc:iec,jsc:jec) )
+allocate( npp3d(isc:iec,jsc:jec,nk) )
 allocate( pprod_gross(isc:iec,jsc:jec,nk) )
 allocate( pprod_gross_2d(isc:iec,jsc:jec) )
 allocate( zprod_gross(isc:iec,jsc:jec,nk) )
@@ -842,7 +909,7 @@ end subroutine  csiro_bgc_end  !}
 !
 
 subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &
-    T_prog, aice, wnd, grid, time, use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, sfc_co2)
+    T_prog, aice, wnd, grid, time, use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, sfc_co2, iof_nit, iof_alg)
 
 use ocmip2_co2calc_mod
 use mpp_mod, only : mpp_sum
@@ -859,6 +926,7 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &
 type(ocean_grid_type), intent(in)                               :: Grid
 type(ocean_time_type), intent(in)                               :: Time
 real, intent(in), dimension(isd:ied,jsd:jed)                    :: aice, wnd
+real, intent(in), dimension(isd:ied,jsd:jed), optional          :: iof_nit, iof_alg
 logical, intent(in) :: use_waterflux, salt_restore_as_salt_flux
 
 real, intent(in), dimension(isd:ied,jsd:jed), optional          :: atm_co2
@@ -947,7 +1015,12 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &
 endif ! if (gasx_from_file)
 
 call time_interp_external(nat_co2_id, time%model_time, nat_co2)
-call time_interp_external(atmpress_id, time%model_time, patm_t)
+if (gasx_from_file) then
+        call time_interp_external(atmpress_id, time%model_time, patm_t)
+else !use the sea level pressure from the forcing (convert Pa to atm)
+        !THIS HAS NOT BEEN IMPLEMENTED YET. SET TO 1 ATM FOR NOW...
+        patm_t(isc:iec,jsc:jec) = 1.0
+endif
 call time_interp_external(dust_id, time%model_time, dust_t)
 if (id_adic .ne. 0) then
 ! The atmospheric co2 value for the anthropogenic+natural carbon tracer
@@ -1271,6 +1344,28 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &
   enddo  !} n
 endif
 
+!ice-to-ocean flux of algae
+if (id_phy.ne.0) then
+  do n = 1, instances  !{
+    do j = jsc, jec  !{
+      do i = isc, iec  !{
+        t_prog(ind_phy)%stf(i,j) = iof_alg(i,j)
+      enddo  !} i
+    enddo  !} j
+  enddo  !} n
+endif
+!ice-to-ocean flux of nitrate
+if (id_no3.ne.0) then
+  do n = 1, instances  !{
+    do j = jsc, jec  !{
+      do i = isc, iec  !{
+        t_prog(ind_no3)%stf(i,j) = iof_nit(i,j)
+      enddo  !} i
+    enddo  !} j
+  enddo  !} n
+endif
+
+
 if (.not. use_waterflux)  then  !{
 ! rjm - One only needs to compute virtual fluxes if waterflux is not used
 ! NB, when this routine is called, t_prog(ind_sal)%stf() only has applied fluxes
@@ -1398,6 +1493,7 @@ subroutine csiro_bgc_init  !{
 real                                                    :: min_range=0.0, max_range=1.e4
 real                                                    :: sum_ntr = 0.0
 character(len=8)                                        :: bgc_trc='tracer00'
+character(len=1)                                        :: bgc_si_prefix='m' ! 'm' for everything but iron, which is 'u'
 
 !       Initialize the csiro_bgc package
 
@@ -1635,12 +1731,19 @@ subroutine csiro_bgc_init  !{
           min_range=0.0
           max_range=100.0
     endif 
+
+    if (nn == id_fe) then
+          bgc_si_prefix = 'u'
+    else
+          bgc_si_prefix='m'
+    endif
           
     
     biotic(n)%ind_bgc(nn) = otpm_set_prog_tracer(trim(bgc_trc) // trim(suffix),     &
        package_name,                                                  &
        longname = trim(bgc_trc) // trim(long_suffix),                      &
-       units = 'mmol/m^3', flux_units = 'mmol/m^2/s',                 &
+       units = bgc_si_prefix//'mol/m^3',                              &
+       flux_units = bgc_si_prefix//'mol/m^2/s',                       &
        caller = trim(mod_name)//'('//trim(sub_name)//')',              &
        min_range=min_range,max_range=max_range)
  enddo  !} nn
@@ -1873,6 +1976,35 @@ subroutine csiro_bgc_source(isc, iec, jsc, jec, isd, ied, jsd, jed, T_prog, grid
        time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
 endif
 
+!det export at 100 m
+if (id_wdet100 .gt. 0) then
+  wdet100(:,:) = wdetbio(isc:iec,jsc:jec)*t_prog(ind_det)%field(isc:iec,jsc:jec,minloc(grid%zt(:)-100,dim=1),time%taum1)
+  used = send_data(id_wdet100, wdet100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+
+! Net primary productivity
+
+! at each depth
+if (id_npp3d .gt. 0) then
+  used = send_data(id_npp3d, npp3d(isc:iec,jsc:jec,:),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,:))
+endif
+! depth integrated
+if (id_npp2d .gt. 0) then
+  npp2d(:,:)=0.0
+  do k=1,grid%nk
+        npp2d(isc:iec,jsc:jec) = npp2d(isc:iec,jsc:jec) + npp3d(isc:iec,jsc:jec,k)*Thickness%dzt(isc:iec,jsc:jec,k)
+  enddo
+  used = send_data(id_npp2d, npp2d(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+! at surface
+if (id_npp1 .gt. 0) then
+  used = send_data(id_npp1, npp3d(isc:iec,jsc:jec,1),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+
 ! Gross production of zooplankton
 
 if (id_zprod_gross .gt. 0) then
@@ -1887,6 +2019,104 @@ subroutine csiro_bgc_source(isc, iec, jsc, jec, isd, ied, jsd, jed, T_prog, grid
        time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
 endif
 
+! mixed-layer-integrated quantities
+
+if (id_adic_intmld .gt. 0) then
+  used = send_data(id_adic_intmld, adic_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_dic_intmld .gt. 0) then
+  used = send_data(id_dic_intmld, dic_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_o2_intmld .gt. 0) then
+  used = send_data(id_o2_intmld, o2_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_no3_intmld .gt. 0) then
+  used = send_data(id_no3_intmld, no3_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_fe_intmld .gt. 0) then
+  used = send_data(id_fe_intmld, fe_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_phy_intmld .gt. 0) then
+  used = send_data(id_phy_intmld, phy_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_det_intmld .gt. 0) then
+  used = send_data(id_det_intmld, det_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_pprod_gross_intmld .gt. 0) then
+  used = send_data(id_pprod_gross_intmld, pprod_gross_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_npp_intmld .gt. 0) then
+  used = send_data(id_npp_intmld, npp_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_radbio_intmld .gt. 0) then
+  used = send_data(id_radbio_intmld, radbio_intmld(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+
+if (id_adic_int100 .gt. 0) then
+  used = send_data(id_adic_int100, adic_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_dic_int100 .gt. 0) then
+  used = send_data(id_dic_int100, dic_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_o2_int100 .gt. 0) then
+  used = send_data(id_o2_int100, o2_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_no3_int100 .gt. 0) then
+  used = send_data(id_no3_int100, no3_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_fe_int100 .gt. 0) then
+  used = send_data(id_fe_int100, fe_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_phy_int100 .gt. 0) then
+  used = send_data(id_phy_int100, phy_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_det_int100 .gt. 0) then
+  used = send_data(id_det_int100, det_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_pprod_gross_int100 .gt. 0) then
+  used = send_data(id_pprod_gross_int100, pprod_gross_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_npp_int100 .gt. 0) then
+  used = send_data(id_npp_int100, npp_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+if (id_radbio_int100 .gt. 0) then
+  used = send_data(id_radbio_int100, radbio_int100(isc:iec,jsc:jec),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+
+! PAR for phytoplankton at surface.
+
+if (id_radbio1 .gt. 0) then
+  used = send_data(id_radbio1, radbio3d(isc:iec,jsc:jec,1),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
+endif
+
+! PAR for phytoplankton at all depths.
+
+if (id_radbio3d .gt. 0) then
+  used = send_data(id_radbio3d, radbio3d(isc:iec,jsc:jec,:),          &
+       time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,:))
+endif
+
 do n = 1, instances  !{
 
  ntr_bgc=biotic(n)%ntr_bgc
@@ -1993,6 +2223,7 @@ subroutine csiro_bgc_start (time, domain, grid) !{
 character(len=5)        :: bgc_btf='btf00'
 character(len=5)        :: bgc_src='src00'
 character(len=64)       :: name1, name2, name3, name4 
+character(len=1)        :: bgc_si_prefix='m' ! 'm' for everything but iron, which is 'u'
 
 
 ! =====================================================================
@@ -2014,23 +2245,27 @@ subroutine csiro_bgc_start (time, domain, grid) !{
 !       Open up the files for boundary conditions
 !-----------------------------------------------------------------------
 
-atmpress_id = init_external_field(atmpress_file,                &
+if (gasx_from_file) then
+   atmpress_id = init_external_field(atmpress_file,                &
                                   atmpress_name,                &
                                   domain = Domain%domain2d)
-if (atmpress_id .eq. 0) then  !{
-  call mpp_error(FATAL, trim(error_header) //                   &
+   if (atmpress_id .eq. 0) then  !{
+      call mpp_error(FATAL, trim(error_header) //                   &
        'Could not open atmpress file: ' //                      &
        trim(atmpress_file))
-endif  !}
+   endif  !}
+endif
 
-pistonveloc_id = init_external_field(pistonveloc_file,          &
+if (gasx_from_file) then
+   pistonveloc_id = init_external_field(pistonveloc_file,          &
                                      pistonveloc_name,          &
                                      domain = Domain%domain2d)
-if (pistonveloc_id .eq. 0) then  !{
-  call mpp_error(FATAL, trim(error_header) //                   &
+   if (pistonveloc_id .eq. 0) then  !{
+     call mpp_error(FATAL, trim(error_header) //                   &
        'Could not open pistonveloc file: ' //                   &
        trim(pistonveloc_file))
-endif  !}
+   endif  !}
+endif
 
 !RASF I think the ifdafs are redundant
 if (id_adic .ne. 0 .and. .not. use_access_co2) then
@@ -2044,14 +2279,16 @@ subroutine csiro_bgc_start (time, domain, grid) !{
  endif  !}
 endif
 
-seaicefract_id = init_external_field(seaicefract_file,          &
+if (ice_file4gasx) then
+   seaicefract_id = init_external_field(seaicefract_file,          &
                                      seaicefract_name,          &
                                      domain = Domain%domain2d)
-if (seaicefract_id .eq. 0) then  !{
-  call mpp_error(FATAL, trim(error_header) //                   &
+   if (seaicefract_id .eq. 0) then  !{
+     call mpp_error(FATAL, trim(error_header) //                   &
        'Could not open seaicefract file: ' //                   &
        trim(seaicefract_file))
-endif  !}
+   endif  !}
+endif
 
 dust_id = init_external_field(dust_file,          &
                                      dust_name,          &
@@ -2203,6 +2440,112 @@ subroutine csiro_bgc_start (time, domain, grid) !{
      grid%tracer_axes(1:2),Time%model_time, 'Integrated light limitation of phytoplankton growth', &
      ' ',missing_value = -1.0e+10)
 
+id_adic_intmld = register_diag_field('ocean_model','adic_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated natural + anthropogenic dissolved inorganic carbon', &
+     'mmol/m^2',missing_value = -1.0e+10)
+id_dic_intmld = register_diag_field('ocean_model','dic_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated natural dissolved inorganic carbon', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_o2_intmld = register_diag_field('ocean_model','o2_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated dissolved oxygen', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_no3_intmld = register_diag_field('ocean_model','no3_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated nitrate', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_fe_intmld = register_diag_field('ocean_model','fe_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated iron', &
+     'umol/m^2',missing_value = -1.0e+10)     
+id_phy_intmld = register_diag_field('ocean_model','phy_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated phytoplankton', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_det_intmld = register_diag_field('ocean_model','det_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated detritus', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_pprod_gross_intmld = register_diag_field('ocean_model','pprod_gross_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated gross PHY production', &
+     'mmolN/m^2/s',missing_value = -1.0e+10)     
+id_npp_intmld = register_diag_field('ocean_model','npp_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated net primary productivity', &
+     'mmolN/m^2/s',missing_value = -1.0e+10)     
+id_radbio_intmld = register_diag_field('ocean_model','radbio_intmld', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     'MLD-integrated photosynthetically active radiation for phytoplankton growth', &
+     'W m-1',missing_value = -1.0e+10)     
+
+id_adic_int100 = register_diag_field('ocean_model','adic_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated natural + anthropogenic dissolved inorganic carbon', &
+     'mmol/m^2',missing_value = -1.0e+10)
+id_dic_int100 = register_diag_field('ocean_model','dic_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated natural dissolved inorganic carbon', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_o2_int100 = register_diag_field('ocean_model','o2_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated dissolved oxygen', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_no3_int100 = register_diag_field('ocean_model','no3_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated nitrate', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_fe_int100 = register_diag_field('ocean_model','fe_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated iron', &
+     'umol/m^2',missing_value = -1.0e+10)     
+id_phy_int100 = register_diag_field('ocean_model','phy_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated phytoplankton', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_det_int100 = register_diag_field('ocean_model','det_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated detritus', &
+     'mmol/m^2',missing_value = -1.0e+10)     
+id_pprod_gross_int100 = register_diag_field('ocean_model','pprod_gross_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated gross PHY production', &
+     'mmolN/m^2/s',missing_value = -1.0e+10)     
+id_npp_int100 = register_diag_field('ocean_model','npp_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated net primary productivity', &
+     'mmolN/m^2/s',missing_value = -1.0e+10)     
+id_radbio_int100 = register_diag_field('ocean_model','radbio_int100', &
+     grid%tracer_axes(1:2),Time%model_time, &
+     '100m-integrated photosynthetically active radiation for phytoplankton growth', &
+     'W m-1',missing_value = -1.0e+10)      
+
+id_radbio1 = register_diag_field('ocean_model','radbio1', &
+     grid%tracer_axes(1:2),Time%model_time, 'Photosynthetically active radiation for phytoplankton growth at surface', &
+     'W m-2',missing_value = -1.0e+10)
+
+id_radbio3d = register_diag_field('ocean_model','radbio3d', &
+     grid%tracer_axes(1:3),Time%model_time, 'Photosynthetically active radiation for phytoplankton growth', &
+     'W m-2',missing_value = -1.0e+10)
+
+id_wdet100 = register_diag_field('ocean_model','wdet100', &
+     grid%tracer_axes(1:2),Time%model_time, 'detritus export at 100 m (det*sinking rate)', &
+     'mmolN/m^2/s',missing_value = -1.0e+10)
+
+id_npp3d = register_diag_field('ocean_model','npp3d', &
+     grid%tracer_axes(1:3),Time%model_time, 'Net primary productivity', &
+     'mmolN/m^3/s',missing_value = -1.0e+10)
+
+id_npp2d = register_diag_field('ocean_model','npp2d', &
+     grid%tracer_axes(1:2),Time%model_time, 'Vertically integrated net primary productivity', &
+     'mmolN/m^2/s',missing_value = -1.0e+10)
+
+id_npp1 = register_diag_field('ocean_model','npp1', &
+     grid%tracer_axes(1:2),Time%model_time, 'Net primary productivity in the first ocean layer', &
+     'mmolN/m^2/s',missing_value = -1.0e+10)
+
 id_pprod_gross = register_diag_field('ocean_model','pprod_gross', &
      grid%tracer_axes(1:3),Time%model_time, 'Gross PHY production', &
      'mmolN/m^3/s',missing_value = -1.0e+10)
@@ -2346,31 +2689,34 @@ subroutine csiro_bgc_start (time, domain, grid) !{
    name2 = 'Virtual flux into ocean - iron'
    name3 = 'Source term - iron'
    name4 = 'Flux into sediment - iron'
+   bgc_si_prefix = 'u'
+  else
+   bgc_si_prefix='m'
   endif
   if (mpp_pe() == mpp_root_pe() )print*,'rjm bio',bgc_stf,'v'//bgc_stf
 
   biotic(n)%id_bgc_stf(nn) = register_diag_field('ocean_model',       &
        bgc_stf//str, grid%tracer_axes(1:2),                     &
-       Time%model_time, name1, 'mmol/m^2/s',    &
-!       Time%model_time, bgc_stf//'flux into ocean', 'mmol/m^2/s',    &
+       Time%model_time, name1, bgc_si_prefix//'mol/m^2/s',    &
+!       Time%model_time, bgc_stf//'flux into ocean', bgc_si_prefix//'mol/m^2/s',    &
        missing_value = -1.0e+10)
 
   biotic(n)%id_bgc_vstf(nn) = register_diag_field('ocean_model',       &
        'v'//bgc_stf//str, grid%tracer_axes(1:2),                     &
-       Time%model_time, name2, 'mmol/m^3/s',    &
-!       Time%model_time, bgc_stf//'virtual flux into ocean', 'mmol/m^3/s',    &
+       Time%model_time, name2, bgc_si_prefix//'mol/m^3/s',    &
+!       Time%model_time, bgc_stf//'virtual flux into ocean', bgc_si_prefix//'mol/m^3/s',    &
        missing_value = -1.0e+10)
 
   biotic(n)%id_bgc_src(nn) = register_diag_field('ocean_model',       &
        bgc_src//str, grid%tracer_axes(1:3),                     &
-       Time%model_time, name3, 'mmolN/m^3/s',    &
-!       Time%model_time, bgc_src, 'mmolN/m^3/s',    &
+       Time%model_time, name3, bgc_si_prefix//'mol/m^3/s',    &
+!       Time%model_time, bgc_src, bgc_si_prefix//'mol/m^3/s',    &
        missing_value = -1.0e+10)
 
   biotic(n)%id_bgc_btf(nn) = register_diag_field('ocean_model',       &
        bgc_btf//str, grid%tracer_axes(1:2),                     &
-       Time%model_time, name4, 'mmol/m^2/s',    &
-!       Time%model_time, bgc_btf//'flux into sediment', 'mmol/m^2/s',    &
+       Time%model_time, name4, bgc_si_prefix//'mol/m^2/s',    &
+!       Time%model_time, bgc_btf//'flux into sediment', bgc_si_prefix//'mol/m^2/s',    &
        missing_value = -1.0e+10)
 
  enddo  !} nn
diff --git a/src/mom5/ocean_tracers/ocean_tpm.F90 b/src/mom5/ocean_tracers/ocean_tpm.F90
index d836419038..576dde39fb 100644
--- a/src/mom5/ocean_tracers/ocean_tpm.F90
+++ b/src/mom5/ocean_tracers/ocean_tpm.F90
@@ -1219,7 +1219,7 @@ end subroutine ocean_tpm_sfc_end  !}
 
 subroutine ocean_tpm_sbc(Domain, Grid, T_prog, Time, Ice_ocean_boundary_fluxes, &
      runoff, isc_bnd, iec_bnd, jsc_bnd, jec_bnd, aice, wnd,            &
-     use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, ocn_co2)
+     use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, ocn_co2, iof_nit, iof_alg)
 
 
 use coupler_types_mod, only: coupler_2d_bc_type
@@ -1247,7 +1247,7 @@ subroutine ocean_tpm_sbc(Domain, Grid, T_prog, Time, Ice_ocean_boundary_fluxes,
 
 real, intent(in), dimension(Domain%isd:,Domain%jsd:), optional :: aice
 real, intent(in), dimension(Domain%isd:,Domain%jsd:), optional :: atm_co2
-real, intent(in), dimension(Domain%isd:,Domain%jsd:), optional :: wnd
+real, intent(in), dimension(Domain%isd:,Domain%jsd:), optional :: wnd, iof_nit, iof_alg
 logical, intent(in), optional                                  :: use_waterflux, salt_restore_as_salt_flux
 
 real, intent(out), dimension(Domain%isd:,Domain%jsd:), optional :: co2flux, ocn_co2
@@ -1321,8 +1321,13 @@ subroutine ocean_tpm_sbc(Domain, Grid, T_prog, Time, Ice_ocean_boundary_fluxes,
 
 #if defined(CSIRO_BGC)
 if (do_csiro_bgc) then  !{
+#if defined(ACCESS_OM)
+  call csiro_bgc_sbc(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, Domain%isd, Domain%ied, Domain%jsd, Domain%jed,  &
+  T_prog, aice, wnd, Grid, Time, use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, ocn_co2, iof_nit=iof_nit, iof_alg=iof_alg)
+#else
   call csiro_bgc_sbc(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, Domain%isd, Domain%ied, Domain%jsd, Domain%jed,  &
   T_prog, aice, wnd, Grid, Time, use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, ocn_co2)
+#endif
 endif  !}
 #endif