Merge pull request mom-ocean#368 from mom-ocean/iamip2-ak

WOMBAT BGC updates

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

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

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)

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 ()