diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90
index a9d4c688f8..1bc4823a20 100644
--- a/config_src/drivers/nuopc_cap/mom_cap.F90
+++ b/config_src/drivers/nuopc_cap/mom_cap.F90
@@ -2,8 +2,9 @@
 
 module MOM_cap_mod
 
+use field_manager_mod,        only: field_manager_init, field_manager_end
 use MOM_domains,              only: get_domain_extent
-use MOM_io,                   only: stdout, io_infra_end
+use MOM_io,                   only: stdout, io_infra_end, slasher
 use mpp_domains_mod,          only: mpp_get_compute_domains
 use mpp_domains_mod,          only: mpp_get_ntile_count, mpp_get_pelist, mpp_get_global_domain
 use mpp_domains_mod,          only: mpp_get_domain_npes
@@ -30,6 +31,7 @@ module MOM_cap_mod
 use MOM_cap_methods,          only: ChkErr
 use MOM_ensemble_manager,     only: ensemble_manager_init
 use MOM_coms,                 only: sum_across_PEs
+use MOM_coupler_types,        only: coupler_1d_bc_type, coupler_2d_bc_type
 
 #ifdef CESMCOUPLED
 use shr_log_mod,             only: shr_log_setLogUnit
@@ -37,6 +39,15 @@ module MOM_cap_mod
 #endif
 use time_utils_mod,           only: esmf2fms_time
 
+#ifdef _USE_GENERIC_TRACER
+use MOM_coupler_types,        only: coupler_type_spawn, coupler_type_destructor
+use MOM_coupler_types,        only: coupler_type_set_diags, coupler_type_send_data, coupler_type_data_override
+use MOM_data_override,        only: data_override_init, data_override
+use MOM_cap_gtracer_flux,     only: gas_exchange_init, gas_fields_restore, gas_fields_restart
+use MOM_cap_gtracer_flux,     only: get_coupled_field_name, add_gas_fluxes_param, UNKNOWN_CMEPS_FIELD
+use MOM_cap_gtracer_flux,     only: atmos_ocean_fluxes_calc
+#endif
+
 use, intrinsic :: iso_fortran_env, only: output_unit
 
 use ESMF,  only: ESMF_ClockAdvance, ESMF_ClockGet, ESMF_ClockPrint, ESMF_VMget
@@ -100,11 +111,14 @@ module MOM_cap_mod
 public SetServices
 public SetVM
 
-!> Internal state type with pointers to three types defined by MOM.
+!> Internal state type with pointers to types defined by MOM.
 type ocean_internalstate_type
   type(ocean_public_type),       pointer :: ocean_public_type_ptr
   type(ocean_state_type),        pointer :: ocean_state_type_ptr
   type(ice_ocean_boundary_type), pointer :: ice_ocean_boundary_type_ptr
+#ifdef _USE_GENERIC_TRACER
+  type(coupler_2d_bc_type),      pointer :: coupler_2d_bc_type_ptr
+#endif
 end type
 
 !>  Wrapper-derived type required to associate an internal state instance
@@ -419,6 +433,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   type (ocean_public_type),      pointer :: ocean_public => NULL()
   type (ocean_state_type),       pointer :: ocean_state => NULL()
   type(ice_ocean_boundary_type), pointer :: Ice_ocean_boundary => NULL()
+  type(coupler_1d_bc_type),      pointer :: gas_fields_atm => NULL()
+  type(coupler_1d_bc_type),      pointer :: gas_fields_ocn => NULL()
+  type(coupler_1d_bc_type),      pointer :: gas_fluxes => NULL()
+  type(coupler_2d_bc_type),      pointer :: atm_fields => NULL()
   type(ocean_internalstate_wrapper)      :: ocean_internalstate
   type(ocean_grid_type),         pointer :: ocean_grid => NULL()
   type(directories)                      :: dirs
@@ -450,6 +468,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   character(len=512)                     :: restartfile          ! Path/Name of restart file
   character(len=2048)                    :: restartfiles         ! Path/Name of restart files
                                                                  ! (same as restartfile if single restart file)
+  character(240)                         :: additional_restart_dir
   character(len=*), parameter            :: subname='(MOM_cap:InitializeAdvertise)'
   character(len=32)                      :: calendar
   character(len=:), allocatable          :: rpointer_filename
@@ -528,6 +547,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   if (chkerr(rc,__LINE__,u_FILE_u)) return
   call MOM_infra_init(mpi_comm_mom)
 
+  call field_manager_init
+
   ! determine the calendar
   if (cesm_coupled) then
     call NUOPC_CompAttributeGet(gcomp, name="calendar", value=cvalue, &
@@ -659,13 +680,44 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
 
   endif
 
+  ! Set NUOPC attribute additional_restart_dir to RESTART/ if not defined
+  additional_restart_dir = "RESTART/"
+  call NUOPC_CompAttributeGet(gcomp, name="additional_restart_dir", value=cvalue, &
+        isPresent=isPresent, isSet=isSet, rc=rc)
+  if (ChkErr(rc,__LINE__,u_FILE_u)) return
+  if (isPresent .and. isSet) then
+    additional_restart_dir = slasher(cvalue)
+  else
+    call ESMF_LogWrite('MOM_cap:additional_restart_dir unset. Defaulting to '//trim(additional_restart_dir), &
+          ESMF_LOGMSG_INFO)
+  endif
+  call NUOPC_CompAttributeSet(gcomp, name="additional_restart_dir", value=additional_restart_dir, rc=rc)
+  if (chkerr(rc,__LINE__,u_FILE_u)) return
+
   ocean_public%is_ocean_pe = .true.
+#ifdef _USE_GENERIC_TRACER
+  ! Initialise structures for extra tracer fluxes
+  call gas_exchange_init(gas_fields_atm=gas_fields_atm, gas_fields_ocn=gas_fields_ocn, gas_fluxes=gas_fluxes)
+
+  if (cesm_coupled .and. len_trim(inst_suffix)>0) then
+    call ocean_model_init(ocean_public, ocean_state, time0, time_start, gas_fields_ocn=gas_fields_ocn, &
+          input_restart_file=trim(adjustl(restartfiles)), inst_index=inst_index)
+  else
+    call ocean_model_init(ocean_public, ocean_state, time0, time_start, gas_fields_ocn=gas_fields_ocn, &
+          input_restart_file=trim(adjustl(restartfiles)))
+  endif
+
+  ! Enable data override via the data_table using the component name 'OCN'
+  call get_ocean_grid(ocean_state, ocean_grid)
+  call data_override_init(ocean_grid%Domain)
+#else
   if (cesm_coupled .and. len_trim(inst_suffix)>0) then
     call ocean_model_init(ocean_public, ocean_state, time0, time_start, &
       input_restart_file=trim(adjustl(restartfiles)), inst_index=inst_index)
   else
     call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(adjustl(restartfiles)))
   endif
+#endif
 
   ! GMM, this call is not needed in CESM. Check with EMC if it can be deleted.
   call ocean_model_flux_init(ocean_state)
@@ -732,6 +784,31 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
     Ice_ocean_boundary%hcond           = 0.0
   endif
 
+#ifdef _USE_GENERIC_TRACER
+  ! Allocate fields for extra tracer fluxes in Ice_ocean_boundary
+  ! Annoyingly, spawning doesn't copy param array, so add manually
+  call coupler_type_spawn(gas_fluxes, Ice_ocean_boundary%fluxes, (/isc,isc,iec,iec/), &
+        (/jsc,jsc,jec,jec/), suffix='_ice_ocn')
+  call add_gas_fluxes_param(Ice_ocean_boundary%fluxes)
+
+  ! Initialise structure for atmos fields related to extra tracer fluxes
+  ! This is set in the ESMF Internal State to be accessed elsewhere
+  ! TODO: should we deallocate atm_fields in a finalise step? Ice_ocean_boundary is handled
+  ! in a similar way and does not appear to be deallocated.
+  allocate(atm_fields)
+  ocean_internalstate%ptr%coupler_2d_bc_type_ptr => atm_fields
+  call coupler_type_spawn(gas_fields_atm, atm_fields, (/isc,isc,iec,iec/), &
+        (/jsc,jsc,jec,jec/), suffix='_atm')
+
+  ! Register diagnosics for extra tracer flux structures
+  call coupler_type_set_diags(Ice_ocean_boundary%fluxes, "ocean_flux", ocean_public%axes(1:2), time_start)
+  call coupler_type_set_diags(atm_fields, "atmos_sfc", ocean_public%axes(1:2), time_start)
+
+  ! Restore ocean fields related to extra tracer fluxes from restart files
+  call get_MOM_input(dirs=dirs)
+  call gas_fields_restore(ocean_public%fields, ocean_public%domain, dirs%restart_input_dir)
+#endif
+
   call query_ocean_state(ocean_state, use_waves=use_waves, wave_method=wave_method)
   if (use_waves) then
     if (wave_method == "EFACTOR") then
@@ -800,6 +877,15 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
     endif
   endif
 
+#ifdef _USE_GENERIC_TRACER
+  ! Add import fields required for extra tracer fluxes
+  do n = 1, gas_fluxes%num_bcs
+    stdname = get_coupled_field_name(gas_fluxes%bc(n)%name)
+    if (stdname /= UNKNOWN_CMEPS_FIELD) &
+      call fld_list_add(fldsToOcn_num, fldsToOcn, stdname, "will provide")
+  enddo
+#endif
+
   !--------- export fields -------------
   call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_omask"   , "will provide")
   call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_t"       , "will provide")
@@ -811,6 +897,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsFrOcn_num, fldsFrOcn, "Fioo_q"     , "will provide")
   call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_bldepth" , "will provide")
 
+  ! TODO: dts: How to handle export fields from generic tracers?
+
   do n = 1,fldsToOcn_num
     call NUOPC_Advertise(importState, standardName=fldsToOcn(n)%stdname, name=fldsToOcn(n)%shortname, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
@@ -1622,11 +1710,14 @@ subroutine ModelAdvance(gcomp, rc)
   type (ocean_public_type),      pointer :: ocean_public       => NULL()
   type (ocean_state_type),       pointer :: ocean_state        => NULL()
   type(ice_ocean_boundary_type), pointer :: Ice_ocean_boundary => NULL()
+  type(coupler_2d_bc_type),      pointer :: atm_fields         => NULL()
   type(ocean_internalstate_wrapper)      :: ocean_internalstate
   type(ocean_grid_type)        , pointer :: ocean_grid
   type(time_type)                        :: Time
+  type(time_type)                        :: Time_import
   type(time_type)                        :: Time_step_coupled
   type(time_type)                        :: Time_restart_current
+  integer                                :: isc,iec,jsc,jec
   integer                                :: dth, dtm, dts
   integer                                :: nc
   type(ESMF_Time)                        :: MyTime
@@ -1638,12 +1729,13 @@ subroutine ModelAdvance(gcomp, rc)
   integer                                :: writeunit
   integer                                :: localPet
   type(ESMF_VM)                          :: vm
-  integer                                :: n, i
+  integer                                :: m, n, i
   character(240)                         :: import_timestr, export_timestr
   character(len=128)                     :: fldname
   character(len=*),parameter             :: subname='(MOM_cap:ModelAdvance)'
   character(len=8)                       :: suffix
   character(len=:), allocatable          :: rpointer_filename
+  character(240)                         :: additional_restart_dir
   integer                                :: num_rest_files
   real(8)                                :: MPI_Wtime, timers
   logical                                :: write_restart, write_restartfh
@@ -1684,6 +1776,7 @@ subroutine ModelAdvance(gcomp, rc)
 
   Time_step_coupled = esmf2fms_time(timeStep)
   Time = esmf2fms_time(currTime)
+  Time_import = Time
 
   !---------------
   ! Apply ocean lag for startup runs:
@@ -1759,8 +1852,39 @@ subroutine ModelAdvance(gcomp, rc)
     ! Import data
     !---------------
 
+#ifdef _USE_GENERIC_TRACER
+    atm_fields => ocean_internalstate%ptr%coupler_2d_bc_type_ptr
+
+    call mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, atm_fields=atm_fields, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    ! Potentially override atm_fields from data_table.
+    call coupler_type_data_override('OCN', atm_fields, Time_import)
+
+    ! Potentially override ice_ocean_boundary%fluxes from data_table.
+    ! Doing this before atmos_ocean_fluxes_calc call avoids unnecessary calculation of overridden fluxes.
+    ! However, we cannot use coupler_type_data_override here since it does not set the override flag on
+    ! overridden fields
+    do n = 1, ice_ocean_boundary%fluxes%num_bcs
+      do m = 1, ice_ocean_boundary%fluxes%bc(n)%num_fields
+        call data_override('OCN', ice_ocean_boundary%fluxes%bc(n)%field(m)%name, &
+              ice_ocean_boundary%fluxes%bc(n)%field(m)%values, Time_import, &
+              override=ice_ocean_boundary%fluxes%bc(n)%field(m)%override)
+      enddo
+    enddo
+
+    ! Calculate the extra tracer fluxes
+    call get_domain_extent(ocean_public%domain, isc, iec, jsc, jec)
+    call atmos_ocean_fluxes_calc(atm_fields, ocean_public%fields, ice_ocean_boundary%fluxes, &
+          ice_ocean_boundary%ice_fraction, isc, iec, jsc, jec)
+
+    ! Send diagnostics
+    call coupler_type_send_data(atm_fields, Time_import)
+    call coupler_type_send_data(ice_ocean_boundary%fluxes, Time_import)
+#else
     call mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
+#endif
 
     !---------------
     ! Update MOM6
@@ -1827,7 +1951,7 @@ subroutine ModelAdvance(gcomp, rc)
       ! determine restart filename
       call ESMF_ClockGetNextTime(clock, MyTime, rc=rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
-      call ESMF_TimeGet (MyTime, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc )
+      call ESMF_TimeGet (MyTime, yy=year, mm=month, dd=day, s=seconds, rc=rc )
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
       if (cesm_coupled) then
@@ -1885,6 +2009,14 @@ subroutine ModelAdvance(gcomp, rc)
 
       endif
 
+#ifdef _USE_GENERIC_TRACER
+      ! Write fields for extra tracer fluxes to their internally defined ocean restart file
+      call NUOPC_CompAttributeGet(gcomp, name="additional_restart_dir", value=additional_restart_dir, rc=rc)
+      if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+      call gas_fields_restart(ocean_public%fields, ocean_public%domain, additional_restart_dir)
+#endif
+
       if (is_root_pe()) then
         write(stdout,*) subname//' writing restart file ',trim(restartname)
       endif
@@ -2124,6 +2256,7 @@ subroutine ocean_model_finalize(gcomp, rc)
   integer                                :: alarmCount
   character(len=64)                      :: timestamp
   logical                                :: write_restart
+  character(240)                         :: additional_restart_dir
   character(len=*),parameter  :: subname='(MOM_cap:ocean_model_finalize)'
   real(8)                                :: MPI_Wtime, timefs
 
@@ -2157,6 +2290,18 @@ subroutine ocean_model_finalize(gcomp, rc)
 
   call ocean_model_end(ocean_public, ocean_State, Time, write_restart=write_restart)
 
+#ifdef _USE_GENERIC_TRACER
+  if (write_restart) then
+    ! Write fields for extra tracer fluxes to their internally defined ocean restart file
+    call NUOPC_CompAttributeGet(gcomp, name="additional_restart_dir", value=additional_restart_dir, rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+
+    call gas_fields_restart(ocean_public%fields, ocean_public%domain, additional_restart_dir)
+  endif
+#endif
+
+  call field_manager_end()
+
   call io_infra_end()
   call MOM_infra_end()
 
diff --git a/config_src/drivers/nuopc_cap/mom_cap_gtracer_flux.F90 b/config_src/drivers/nuopc_cap/mom_cap_gtracer_flux.F90
new file mode 100644
index 0000000000..af1d886b37
--- /dev/null
+++ b/config_src/drivers/nuopc_cap/mom_cap_gtracer_flux.F90
@@ -0,0 +1,680 @@
+!> Contains routines for handling FMS coupler_bc_type tracer flux structures
+!! when using MOM generic tracers
+
+module MOM_cap_gtracer_flux
+
+use MOM_domains,               only: domain2d
+use MOM_coupler_types,         only: coupler_1d_bc_type, coupler_2d_bc_type
+use MOM_coupler_types,         only: ind_flux, ind_deltap, ind_kw, ind_flux0
+use MOM_coupler_types,         only: ind_pcair, ind_u10, ind_psurf
+use MOM_coupler_types,         only: ind_alpha, ind_csurf, ind_sc_no
+use MOM_coupler_types,         only: ind_runoff, ind_deposition
+use MOM_ocean_model_nuopc,     only: ocean_model_flux_init
+use coupler_types_mod,         only: coupler_type_register_restarts, coupler_type_restore_state
+use atmos_ocean_fluxes_mod,    only: atmos_ocean_type_fluxes_init, atmos_ocean_fluxes_init
+use fms2_io_mod,               only: FmsNetcdfDomainFile_t
+use fms2_io_mod,               only: fms2_check_if_open => check_if_open, fms2_close_file => close_file
+use fms2_io_mod,               only: fms2_write_data => write_data
+use fms2_io_mod,               only: fms2_read_restart => read_restart, fms2_write_restart => write_restart
+use fms2_io_mod,               only: fms2_get_global_io_domain_indices => get_global_io_domain_indices
+use field_manager_mod,         only: fm_field_name_len, fm_type_name_len, fm_loop_over_list, fm_change_list
+use fm_util_mod,               only: fm_util_get_real_array
+use mpp_mod,                   only: mpp_error, FATAL
+use FMSconstants,              only: wtmair, rdgas, vonkarm
+
+implicit none; private
+
+! Public member functions
+public :: gas_exchange_init
+public :: gas_fields_restore
+public :: gas_fields_restart
+public :: add_gas_fluxes_param
+public :: get_coupled_field_name
+public :: atmos_ocean_fluxes_calc
+public :: UNKNOWN_CMEPS_FIELD
+
+character(len=*), parameter      :: mod_name = 'mom_cap_gtracer_flux'
+character(len=*), parameter      :: UNKNOWN_CMEPS_FIELD = "UNKNOWN_FIELD"
+real, parameter                  :: epsln=1.0e-30
+
+!> FMS coupler_bc_types for additional tracer fields when using generic tracers
+logical :: gas_fluxes_initialized = .false.  ! This is set to true when the following types are initialized.
+type(coupler_1d_bc_type), target :: ex_gas_fields_atm ! tracer fields in atm
+    !< Structure containing atmospheric surface variables that are used in the
+    !! calculation of the atmosphere-ocean gas fluxes, as well as parameters
+    !! regulating these fluxes. The fields in this structure are never actually
+    !! set, but the structure is used for initialisation of components and to
+    !! spawn other structure whose fields are set.
+type(coupler_1d_bc_type), target :: ex_gas_fields_ocn ! tracer fields atop the ocean
+    !< Structure containing ocean surface variables that are used in the
+    !! calculation of the atmosphere-ocean gas fluxes, as well as parameters
+    !! regulating these fluxes. The fields in this structure are never actually
+    !! set, but the structure is used for initialisation of components and to
+    !! spawn other structure whose fields are set.
+type(coupler_1d_bc_type), target :: ex_gas_fluxes ! tracer fluxes between the atm and ocean
+    !< A structure for exchanging gas or tracer fluxes between the atmosphere and
+    !! ocean, defined by the field table, as well as a place holder of
+    !! intermediate calculations, such as piston velocities, and parameters that
+    !! impact the fluxes. The fields in this structure are never actually set,
+    !! but the structure is used for initialisation of components and to spawn
+    !! other structure whose fields are set.
+
+contains
+
+!> \brief Gas and tracer field initialization routine for running with MOM generic tracers.
+!! Copied and adapted slightly from
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/7761886/full/flux_exchange.F90#L626.
+subroutine gas_exchange_init (gas_fields_atm, gas_fields_ocn, gas_fluxes)
+  type(coupler_1d_bc_type), optional, pointer :: gas_fields_atm ! tracer fields in atm
+      !< Pointer to a structure containing atmospheric surface variables that
+      !! are used in the calculation of the atmosphere-ocean gas fluxes, as well
+      !! as parameters regulating these fluxes.
+  type(coupler_1d_bc_type), optional, pointer :: gas_fields_ocn ! tracer fields atop the ocean
+      !< Pointer to a structure containing ocean surface variables that are
+      !! used in the calculation of the atmosphere-ocean gas fluxes, as well as
+      !! parameters regulating these fluxes.
+  type(coupler_1d_bc_type), optional, pointer :: gas_fluxes ! tracer fluxes between the atm and ocean
+      !< Pointer to a structure for exchanging gas or tracer fluxes between the
+      !! atmosphere and ocean, defined by the field table, as well as a place holder
+      !! of intermediate calculations, such as piston velocities, and parameters
+      !! that impact the fluxes.
+
+  if (.not.gas_fluxes_initialized) then
+    call atmos_ocean_type_fluxes_init( )
+    call ocean_model_flux_init( )
+    call atmos_ocean_fluxes_init(ex_gas_fluxes, ex_gas_fields_atm, ex_gas_fields_ocn)
+    gas_fluxes_initialized = .true.
+  endif
+
+  if (present(gas_fields_atm)) gas_fields_atm => ex_gas_fields_atm
+  if (present(gas_fields_ocn)) gas_fields_ocn => ex_gas_fields_ocn
+  if (present(gas_fluxes)) gas_fluxes => ex_gas_fluxes
+
+end subroutine gas_exchange_init
+
+!> \brief Restore FMS coupler_bc_type state from ocean restart file
+! See https://github.com/NOAA-GFDL/FMScoupler/blob/008399d/full/full_coupler_mod.F90#L1096
+subroutine gas_fields_restore(gas_fields, domain, directory)
+  type(coupler_2d_bc_type), intent(inout) :: gas_fields !< FMS coupler_bc_type to be registered for restarts
+  type(domain2D), intent(in)              :: domain     !< The FMS domain to use for this registration call
+  character(len=*), optional, intent(in)  :: directory  !< Directory containing the restart file
+
+  ! local variables
+  type(FmsNetcdfDomainFile_t), dimension(:), pointer :: ocn_bc_restart => NULL() !< Structures describing
+                                                                                 !! the restart files
+  integer                                            :: num_ocn_bc_restart       !< The number of restart
+                                                                                 !! files to use
+  integer                                            :: n
+
+  call coupler_type_register_restarts(gas_fields, ocn_bc_restart, num_ocn_bc_restart, &
+      domain, to_read=.true., ocean_restart=.true., directory=directory)
+
+  ! Restore the fields from the restart files
+  do n = 1, num_ocn_bc_restart
+    if (fms2_check_if_open(ocn_bc_restart(n))) then
+      call fms2_read_restart(ocn_bc_restart(n))
+    endif
+  enddo
+
+  ! Check whether the restarts were read successfully.
+  call coupler_type_restore_state(gas_fields, use_fms2_io=.true., test_by_field=.true.)
+
+  do n = 1, num_ocn_bc_restart
+    if(fms2_check_if_open(ocn_bc_restart(n))) call fms2_close_file(ocn_bc_restart(n))
+  enddo
+
+end subroutine gas_fields_restore
+
+!> \brief Write ocean restart file for FMS coupler_bc_type state
+! See https://github.com/NOAA-GFDL/FMScoupler/blob/008399d/full/full_coupler_mod.F90#L1408
+subroutine gas_fields_restart(gas_fields, domain, directory)
+  type(coupler_2d_bc_type), intent(inout) :: gas_fields !< FMS coupler_bc_type to be registered for restarts
+  type(domain2D), intent(in)              :: domain     !< The FMS domain to use for this registration call
+  character(len=*), optional, intent(in)  :: directory  !< Directory containing the restart file
+
+  ! local variables
+  type(FmsNetcdfDomainFile_t), dimension(:), pointer :: ocn_bc_restart => NULL() !< Structures describing
+                                                                                 !! the restart files
+  integer                                            :: num_ocn_bc_restart       !< The number of restart
+                                                                                 !! files to use
+  integer                                            :: n
+
+  call coupler_type_register_restarts(gas_fields, ocn_bc_restart, num_ocn_bc_restart, &
+      domain, to_read=.false., ocean_restart=.true., directory=directory)
+
+  do n = 1, num_ocn_bc_restart
+    if (fms2_check_if_open(ocn_bc_restart(n))) then
+      call fms2_write_restart(ocn_bc_restart(n))
+      call add_domain_dimension_data(ocn_bc_restart(n))
+      call fms2_close_file(ocn_bc_restart(n))
+    endif
+  enddo
+
+end subroutine gas_fields_restart
+
+!> Register the axis data as a variable in the netcdf file and add some dummy data.
+!! This is needed so the combiner can work correctly when the io_layout is not 1,1. Copied from
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/008399d/full/full_coupler_mod.F90#L1328
+subroutine add_domain_dimension_data(fileobj)
+  type(FmsNetcdfDomainFile_t)        :: fileobj !< Fms2io domain decomposed fileobj
+
+  ! local variables
+  integer, dimension(:), allocatable :: buffer !< Buffer with axis data
+  integer                            :: is, ie !< Starting and Ending indices for data
+
+  call fms2_get_global_io_domain_indices(fileobj, "xaxis_1", is, ie, indices=buffer)
+  call fms2_write_data(fileobj, "xaxis_1", buffer)
+  deallocate(buffer)
+
+  call fms2_get_global_io_domain_indices(fileobj, "yaxis_1", is, ie, indices=buffer)
+  call fms2_write_data(fileobj, "yaxis_1", buffer)
+  deallocate(buffer)
+
+end subroutine add_domain_dimension_data
+
+!> Retrieve param array from field_manager and add to FMS coupler_bc_type. This is
+!! needed because the coupler_type_spawn routine does not copy the param array into
+!! the spawned type. Hopefully we can get rid of this. This routine is based on
+!! https://github.com/NOAA-GFDL/FMS/blob/7f58528/coupler/atmos_ocean_fluxes.F90#L448
+subroutine add_gas_fluxes_param(gas_fluxes)
+  type(coupler_2d_bc_type), intent(inout) :: gas_fluxes !< FMS coupler_bc_type to add param to
+
+  ! local variables
+  integer                                 :: n
+  character(len=fm_field_name_len)        :: name
+  character(len=fm_type_name_len)         :: typ
+  integer                                 :: ind
+
+  character(len=*), parameter             :: sub_name = 'add_gas_fluxes_param'
+  character(len=*), parameter             :: error_header =&
+      '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
+
+  n = 0
+  do while (fm_loop_over_list('/coupler_mod/fluxes', name, typ, ind))
+    if (typ .ne. 'list') then
+      call mpp_error(FATAL, trim(error_header) // ' ' // trim(name) // ' is not a list')
+    endif
+
+    n = n + 1
+
+    if (.not. fm_change_list('/coupler_mod/fluxes/' // trim(name))) then
+      call mpp_error(FATAL, trim(error_header) // ' Problem changing to ' // trim(name))
+    endif
+
+    if (gas_fluxes%bc(n)%name .eq. name) then
+      gas_fluxes%bc(n)%param => fm_util_get_real_array('param')
+    else
+      call mpp_error(FATAL, trim(error_header) // ' Problem setting param array pointer')
+    endif
+  enddo
+end subroutine add_gas_fluxes_param
+
+!> Return the CMEPS standard_name of the coupled field required for a given coupled
+!! generic_tracer flux name.
+function get_coupled_field_name(name)
+  character(len=64)                  :: get_coupled_field_name !< CMEPS standard_name
+  character(len=*), intent(in)       :: name                   !< gtracer flux name
+
+  ! Add other coupled field names here
+  select case(trim(name))
+    case( 'co2_flux' )
+      get_coupled_field_name = "Sa_co2prog"
+    case default
+      get_coupled_field_name = UNKNOWN_CMEPS_FIELD
+  end select
+end function get_coupled_field_name
+
+!> \brief Calculate the FMS coupler_bc_type ocean tracer fluxes. Units should be mol/m^2/s.
+!! Upward flux is positive.
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+!! and subsequently modified in the following ways:
+!! - Operate on 2D inputs, rather than 1D
+!! - Add calculation for 'air_sea_deposition' taken from
+!!   https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_dep_fluxes_calc.F90
+!! - Multiply fluxes by ice_fraction input, rather than masking based on seawater input
+!! - Use MOM over FMS modules where easy to do so
+!! - Make tsurf input optional, as it is only used by a few implementations
+!! - Use ind_runoff rather than ind_deposition in runoff flux calculation (note, their
+!!   values are equal)
+!! - Rename gas_fields_ice to gas_fields_ocn
+subroutine atmos_ocean_fluxes_calc(gas_fields_atm, gas_fields_ocn, gas_fluxes,&
+  ice_fraction, isc, iec, jsc, jec, tsurf, ustar, cd_m)
+  type(coupler_2d_bc_type), intent(in)     :: gas_fields_atm ! fields in atm
+      !< Structure containing atmospheric surface variables that are used in the calculation
+      !! of the atmosphere-ocean tracer fluxes.
+  type(coupler_2d_bc_type), intent(in)     :: gas_fields_ocn ! fields atop the ocean
+      !< Structure containing ocean surface variables that are used in the calculation of the
+      !! atmosphere-ocean tracer fluxes.
+  type(coupler_2d_bc_type), intent(inout)  :: gas_fluxes ! fluxes between the atm and ocean
+      !< Structure containing the gas fluxes between the atmosphere and the ocean and
+      !! parameters related to the calculation of these fluxes.
+  real, intent(in)                         :: ice_fraction(isc:iec,jsc:jec) !< sea ice fraction
+  integer, intent(in)                      :: isc !< The start i-index of cell centers within
+                                                  !! the computational domain
+  integer, intent(in)                      :: iec !< The end i-index of cell centers within the
+                                                  !! computational domain
+  integer, intent(in)                      :: jsc !< The start j-index of cell centers within
+                                                  !! the computational domain
+  integer, intent(in)                      :: jec !< The end j-index of cell centers within the
+                                                  !! computational domain
+  real, intent(in), optional               :: tsurf(isc:iec,jsc:jec) !< surface temperature
+  real, intent(in), optional               :: ustar(isc:iec,jsc:jec) !< friction velocity, not
+                                                                      !! used
+  real, intent(in), optional               :: cd_m (isc:iec,jsc:jec) !< drag coefficient, not
+                                                                      !! used
+
+  ! local variables
+  character(len=*), parameter   :: sub_name = 'atmos_ocean_fluxes_calc'
+  character(len=*), parameter   :: error_header =&
+      & '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
+  real, parameter                         :: permeg=1.0e-6
+
+  integer                                 :: n
+  integer                                 :: i
+  integer                                 :: j
+  real, dimension(:,:), allocatable       :: kw
+  real, dimension(:,:), allocatable       :: cair
+  character(len=128)                      :: error_string
+
+  ! Return if no fluxes to be calculated
+  if (gas_fluxes%num_bcs .le. 0) return
+
+  if (.not. associated(gas_fluxes%bc)) then
+    if (gas_fluxes%num_bcs .ne. 0) then
+      call mpp_error(FATAL, trim(error_header) // ' Number of gas fluxes not zero')
+    else
+      return
+    endif
+  endif
+
+  do n = 1, gas_fluxes%num_bcs
+    ! only do calculations if the flux has not been overridden
+    if ( .not. gas_fluxes%bc(n)%field(ind_flux)%override) then
+      if (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_gas_flux_generic') then
+        if (.not. allocated(kw)) then
+          allocate( kw(isc:iec,jsc:jec) )
+          allocate ( cair(isc:iec,jsc:jec) )
+        elseif ((size(kw(:,:), dim=1) .ne. iec-isc+1) .or. (size(kw(:,:), dim=2) .ne. jec-jsc+1)) then
+          call mpp_error(FATAL, trim(error_header) // ' Sizes of flux fields do not match')
+        endif
+
+        if (gas_fluxes%bc(n)%implementation .eq. 'ocmip2') then
+          do j = jsc,jec
+            do i = isc,iec
+              gas_fluxes%bc(n)%field(ind_kw)%values(i,j) =&
+                  & (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) * &
+                  & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j)**2
+              cair(i,j) = &
+                  gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * &
+                  gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) * &
+                  gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2)
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =&
+                  & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *&
+                  & sqrt(660. / (gas_fields_ocn%bc(n)%field(ind_sc_no)%values(i,j) + epsln)) *&
+                  & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j))
+              gas_fluxes%bc(n)%field(ind_flux0)%values(i,j) =&
+                  & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *&
+                  & sqrt(660. / (gas_fields_ocn%bc(n)%field(ind_sc_no)%values(i,j) + epsln)) *&
+                  & gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j)
+              gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =&
+                  & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) / &
+                (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln)
+            enddo
+          enddo
+        elseif (gas_fluxes%bc(n)%implementation .eq. 'duce') then
+          if (.not. present(tsurf)) then
+            call mpp_error(FATAL, trim(error_header) // ' Implementation ' //&
+                trim(gas_fluxes%bc(n)%implementation) // ' for ' // trim(gas_fluxes%bc(n)%name) //&
+                ' requires input tsurf')
+          endif
+          do j = jsc,jec
+            do i = isc,iec
+              gas_fluxes%bc(n)%field(ind_kw)%values(i,j) = &
+                  & (1 - ice_fraction(i,j)) * gas_fields_atm%bc(n)%field(ind_u10)%values(i,j) /&
+                  & (770.+45.*gas_fluxes%bc(n)%param(1)**(1./3.)) *&
+                  & 101325./(rdgas*wtmair*1e-3*tsurf(i,j) *&
+                  & max(gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j),epsln))
+              !alpha: mol/m3/atm
+              cair(i,j) = &
+                  gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * &
+                  gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) * &
+                  gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * 9.86923e-6
+              cair(i,j) = max(cair(i,j),0.)
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =&
+                  & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *&
+                  & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j))
+              gas_fluxes%bc(n)%field(ind_flux0)%values(i,j) =&
+                  & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *&
+                  & max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.)
+              gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =&
+                  & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) /&
+                  & (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln)
+            enddo
+          enddo
+        elseif (gas_fluxes%bc(n)%implementation .eq. 'johnson') then
+          if (.not. present(tsurf)) then
+            call mpp_error(FATAL, trim(error_header) // ' Implementation ' //&
+                trim(gas_fluxes%bc(n)%implementation) // ' for ' // trim(gas_fluxes%bc(n)%name) //&
+                ' requires input tsurf')
+          endif
+          !f1p: not sure how to pass salinity. For now, just force at 35.
+          do j = jsc,jec
+            do i = isc,iec
+              !calc_kw(tk,p,u10,h,vb,mw,sc_w,ustar,cd_m)
+              gas_fluxes%bc(n)%field(ind_kw)%values(i,j) =&
+                  & (1 - ice_fraction(i,j)) * calc_kw(tsurf(i,j),&
+                  & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j),&
+                  & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j),&
+                  & 101325./(rdgas*wtmair*1e-3*tsurf(i,j)*&
+                  & max(gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j),epsln)),&
+                  & gas_fluxes%bc(n)%param(2),&
+                  & gas_fluxes%bc(n)%param(1),&
+                  & gas_fields_ocn%bc(n)%field(ind_sc_no)%values(i,j))
+              cair(i,j) =&
+                  & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *&
+                  & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *&
+                  & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * 9.86923e-6
+              cair(i,j) = max(cair(i,j),0.)
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =&
+                  & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *&
+                  & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j))
+              gas_fluxes%bc(n)%field(ind_flux0)%values(i,j) =&
+                  & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *&
+                  & max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.)
+              gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =&
+                  & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) /&
+                  & (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln)
+            enddo
+          enddo
+        else
+          call mpp_error(FATAL, ' Unknown implementation (' //&
+              & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name))
+        endif
+      elseif (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_gas_flux') then
+        if (.not. allocated(kw)) then
+          allocate( kw(isc:iec,jsc:jec) )
+          allocate ( cair(isc:iec,jsc:jec) )
+        elseif ((size(kw(:,:), dim=1) .ne. iec-isc+1) .or. (size(kw(:,:), dim=2) .ne. jec-jsc+1)) then
+          call mpp_error(FATAL, trim(error_header) // ' Sizes of flux fields do not match')
+        endif
+
+        if (gas_fluxes%bc(n)%implementation .eq. 'ocmip2_data') then
+          do j = jsc,jec
+            do i = isc,iec
+              kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *&
+                  & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j)
+              cair(i,j) =&
+                  & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *&
+                  & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *&
+                  & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2)
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *&
+                  & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j))
+            enddo
+          enddo
+        elseif (gas_fluxes%bc(n)%implementation .eq. 'ocmip2') then
+          do j = jsc,jec
+            do i = isc,iec
+              kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *&
+                  & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j)**2
+              cair(i,j) =&
+                  & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *&
+                  & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *&
+                  & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2)
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *&
+                  & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j))
+            enddo
+          enddo
+        elseif (gas_fluxes%bc(n)%implementation .eq. 'linear') then
+          do j = jsc,jec
+            do i = isc,iec
+              kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *&
+                  & max(0.0, gas_fields_atm%bc(n)%field(ind_u10)%values(i,j) - gas_fluxes%bc(n)%param(2))
+              cair(i,j) =&
+                  & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *&
+                  & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *&
+                  & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(3)
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *&
+                  & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j))
+            enddo
+          enddo
+        else
+          call mpp_error(FATAL, ' Unknown implementation (' //&
+              & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name))
+        endif
+      elseif (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_deposition') then
+        if (gas_fluxes%bc(n)%param(1) .le. 0.0) then
+          write (error_string, '(1pe10.3)') gas_fluxes%bc(n)%param(1)
+          call mpp_error(FATAL, 'Bad parameter (' // trim(error_string) //&
+              & ') for air_sea_deposition for ' // trim(gas_fluxes%bc(n)%name))
+        endif
+
+        if (gas_fluxes%bc(n)%implementation .eq. 'dry') then
+          do j = jsc,jec
+            do i = isc,iec
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *&
+                  gas_fields_atm%bc(n)%field(ind_deposition)%values(i,j) / gas_fluxes%bc(n)%param(1)
+            enddo
+          enddo
+        elseif (gas_fluxes%bc(n)%implementation .eq. 'wet') then
+          do j = jsc,jec
+            do i = isc,iec
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *&
+                  gas_fields_atm%bc(n)%field(ind_deposition)%values(i,j) / gas_fluxes%bc(n)%param(1)
+            enddo
+          enddo
+        else
+          call mpp_error(FATAL, 'Unknown implementation (' //&
+              & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name))
+        endif
+      elseif (gas_fluxes%bc(n)%flux_type .eq. 'land_sea_runoff') then
+        if (gas_fluxes%bc(n)%param(1) .le. 0.0) then
+          write (error_string, '(1pe10.3)') gas_fluxes%bc(n)%param(1)
+          call mpp_error(FATAL, ' Bad parameter (' // trim(error_string) //&
+              & ') for land_sea_runoff for ' // trim(gas_fluxes%bc(n)%name))
+        endif
+
+        if (gas_fluxes%bc(n)%implementation .eq. 'river') then
+          do j = jsc,jec
+            do i = isc,iec
+              gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *&
+                  & gas_fields_atm%bc(n)%field(ind_runoff)%values(i,j) /&
+                  & gas_fluxes%bc(n)%param(1)
+            enddo
+          enddo
+        else
+          call mpp_error(FATAL, ' Unknown implementation (' //&
+              & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name))
+        endif
+      else
+        call mpp_error(FATAL, ' Unknown flux_type (' // trim(gas_fluxes%bc(n)%flux_type) //&
+            & ') for ' // trim(gas_fluxes%bc(n)%name))
+      endif
+    endif
+  enddo
+
+  if (allocated(kw)) then
+    deallocate(kw)
+    deallocate(cair)
+  endif
+end subroutine  atmos_ocean_fluxes_calc
+
+!> Calculate \f$k_w\f$
+!!
+!! Taken from Johnson, Ocean Science, 2010. (http://doi.org/10.5194/os-6-913-2010)
+!!
+!! Uses equations defined in Liss[1974],
+!! \f[
+!!  F = K_g(c_g - H C_l) = K_l(c_g/H - C_l)
+!! \f]
+!! where \f$c_g\f$ and \f$C_l\f$ are the bulk gas and liquid concentrations, \f$H\f$
+!! is the Henry's law constant (\f$H = c_{sg}/C_{sl}\f$, where \f$c_{sg}\f$ is the
+!! equilibrium concentration in gas phase (\f$g/cm^3\f$ of air) and \f$C_{sl}\f$ is the
+!! equilibrium concentration of unionised dissolved gas in liquid phase (\f$g/cm^3\f$
+!! of water)),
+!! \f[
+!!    1/K_g = 1/k_g + H/k_l
+!! \f]
+!! and
+!! \f[
+!!    1/K_l = 1/k_l + 1/{Hk_g}
+!! \f]
+!! where \f$k_g\f$ and \f$k_l\f$ are the exchange constants for the gas and liquid
+!! phases, respectively.
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+real function calc_kw(tk, p, u10, h, vb, mw, sc_w, ustar, cd_m)
+  real, intent(in) :: tk !< temperature at surface in kelvin
+  real, intent(in) :: p !< pressure at surface in pa
+  real, intent(in) :: u10 !< wind speed at 10m above the surface in m/s
+  real, intent(in) :: h !< Henry's law constant (\f$H=c_sg/C_sl\f$) (unitless)
+  real, intent(in) :: vb !< Molar volume
+  real, intent(in) :: mw !< molecular weight (g/mol)
+  real, intent(in) :: sc_w
+  real, intent(in), optional :: ustar !< Friction velocity (m/s).  If not provided,
+                                      !! ustar = \f$u_{10} \sqrt{C_D}\f$.
+  real, intent(in), optional :: cd_m !< Drag coefficient (\f$C_D\f$).  Used only if
+                                      !! ustar is provided.
+                                      !! If ustar is not provided,
+                                      !! cd_m = \f$6.1 \times 10^{-4} + 0.63 \times 10^{-4} *u_10\f$
+
+  real :: ra,rl,tc
+
+  tc = tk-273.15
+  ra = 1./max(h*calc_ka(tc,p,mw,vb,u10,ustar,cd_m),epsln)
+  rl = 1./max(calc_kl(tc,u10,sc_w),epsln)
+  calc_kw = 1./max(ra+rl,epsln)
+end function calc_kw
+
+!> Calculate \f$k_a\f$
+!!
+!! See calc_kw
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+real function calc_ka(t, p, mw, vb, u10, ustar, cd_m)
+  real, intent(in) :: t !< temperature at surface in C
+  real, intent(in) :: p !< pressure at surface in pa
+  real, intent(in) :: mw !< molecular weight (g/mol)
+  real, intent(in) :: vb !< molar volume
+  real, intent(in) :: u10 !< wind speed at 10m above the surface in m/s
+  real, intent(in), optional :: ustar !< Friction velocity (m/s).  If not provided,
+                                      !! ustar = \f$u_{10} \sqrt{C_D}\f$.
+  real, intent(in), optional :: cd_m !< Drag coefficient (\f$C_D\f$).  Used only if
+                                      !! ustar is provided.
+                                      !! If ustar is not provided,
+                                      !! cd_m = \f$6.1 \times 10^{-4} + 0.63 \times 10^{-4} *u_10\f$
+
+  real             :: sc
+  real             :: ustar_t, cd_m_t
+
+  if (.not. present(ustar)) then
+    !drag coefficient
+    cd_m_t = 6.1e-4 +0.63e-4*u10
+    !friction velocity
+    ustar_t = u10*sqrt(cd_m_t)
+  else
+    cd_m_t = cd_m
+    ustar_t = ustar
+  end if
+  sc = schmidt_g(t,p,mw,vb)
+  calc_ka = 1e-3+ustar_t/(13.3*sqrt(sc)+1/sqrt(cd_m_t)-5.+log(sc)/(2.*vonkarm))
+end function calc_ka
+
+!> Calculate \f$k_l\f$
+!!
+!! See calc_kw, and Nightingale, Global Biogeochemical Cycles, 2000
+!! (https://doi.org/10.1029/1999GB900091)
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+real function calc_kl(t, v, sc)
+  real, intent(in) :: t !< temperature at surface in C
+  real, intent(in) :: v !< wind speed at surface in m/s
+  real, intent(in) :: sc
+
+  calc_kl = (((0.222*v**2)+0.333*v)*(max(sc,epsln)/600.)**(-0.5))/(100.*3600.)
+end function calc_kl
+
+!> Schmidt number of the gas in air
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+real function schmidt_g(t, p, mw, vb)
+  real, intent(in) :: t !< temperature at surface in C
+  real, intent(in) :: p !< pressure at surface in pa
+  real, intent(in) :: mw !< molecular weight (g/mol)
+  real, intent(in) :: vb !< molar volume
+
+  real :: d,v
+
+  d = d_air(t,p,mw,vb)
+  v = v_air(t)
+  schmidt_g = v / d
+end function schmidt_g
+
+!> From Fuller, Industrial & Engineering Chemistry (https://doi.org/10.1021/ie50677a007)
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+real function d_air(t, p, mw, vb)
+  real, intent(in) :: t  !< temperature in c
+  real, intent(in) :: p  !< pressure in pa
+  real, intent(in) :: mw !< molecular weight (g/mol)
+  real, intent(in) :: vb !< diffusion coefficient (\f$cm3/mol\f$)
+
+  real, parameter :: ma = 28.97d0 !< molecular weight air in g/mol
+  real, parameter :: va = 20.1d0  !< diffusion volume for air (\f$cm^3/mol\f$)
+
+  real            :: pa
+
+  ! convert p to atm
+  pa = 9.8692d-6*p
+  d_air = 1d-3 *&
+      & (t+273.15d0)**(1.75d0)*sqrt(1d0/ma + 1d0/mw)/(pa*(va**(1d0/3d0)+vb**(1d0/3d0))**2d0)
+  ! d_air is in cm2/s convert to m2/s
+  d_air = d_air * 1d-4
+end function d_air
+
+!> kinematic viscosity in air
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+real function p_air(t)
+  real, intent(in) :: t
+
+  real, parameter :: sd_0 = 1.293393662d0,&
+      & sd_1 = -5.538444326d-3,&
+      & sd_2 = 3.860201577d-5,&
+      & sd_3 = -5.2536065d-7
+  p_air = sd_0+(sd_1*t)+(sd_2*t**2)+(sd_3*t**3)
+end function p_air
+
+!> Kinematic viscosity in air (\f$m^2/s\f$
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+real function v_air(t)
+  real, intent(in) :: t !< temperature in C
+  v_air = n_air(t)/p_air(t)
+end function v_air
+
+!> dynamic viscosity in air
+!!
+!! This routine was copied from FMScoupler at
+!! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d38/full/atmos_ocean_fluxes_calc.F90
+real function n_air(t)
+  real, intent(in) :: t !< temperature in C
+
+  real, parameter :: sv_0 = 1.715747771d-5,&
+      & sv_1 = 4.722402075d-8,&
+      & sv_2 = -3.663027156d-10,&
+      & sv_3 = 1.873236686d-12,&
+      & sv_4 = -8.050218737d-14
+  ! in n.s/m^2 (pa.s)
+  n_air = sv_0+(sv_1*t)+(sv_2*t**2)+(sv_3*t**3)+(sv_4*t**4)
+end function n_air
+
+end module MOM_cap_gtracer_flux
\ No newline at end of file
diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90
index 125bae5748..3f81f6aa7e 100644
--- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90
+++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90
@@ -20,8 +20,15 @@ module MOM_cap_methods
 use MOM_surface_forcing_nuopc, only: ice_ocean_boundary_type
 use MOM_grid,                  only: ocean_grid_type
 use MOM_domains,               only: pass_var
+use MOM_coupler_types,         only: coupler_2d_bc_type
 use mpp_domains_mod,           only: mpp_get_compute_domain
 
+#ifdef _USE_GENERIC_TRACER
+use MOM_coupler_types,         only: set_coupler_type_data
+use MOM_coupler_types,         only: ind_pcair, ind_u10, ind_psurf, ind_runoff, ind_deposition
+use MOM_cap_gtracer_flux,      only: get_coupled_field_name, UNKNOWN_CMEPS_FIELD
+#endif
+
 ! By default make data private
 implicit none; private
 
@@ -72,11 +79,17 @@ end subroutine mom_set_geomtype
 !> This function has a few purposes:
 !! (1) it imports surface fluxes using data from the mediator; and
 !! (2) it can apply restoring in SST and SSS.
-subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, rc)
+!! (3) optional: if atm_fields is provided, it imports and sets the fields in atm_fields required
+!! for the calculation of coupled generic tracer fluxes
+subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, atm_fields, rc)
   type(ocean_public_type)       , intent(in)    :: ocean_public       !< Ocean surface state
   type(ocean_grid_type)         , intent(in)    :: ocean_grid         !< Ocean model grid
   type(ESMF_State)              , intent(inout) :: importState        !< incoming data from mediator
   type(ice_ocean_boundary_type) , intent(inout) :: ice_ocean_boundary !< Ocean boundary forcing
+  type(coupler_2d_bc_type), optional, intent(inout) :: atm_fields     !< If present, this type
+                                                       !! describes the atmospheric tracer fields to
+                                                       !! be imported for the calculation of generic
+                                                       !! tracer fluxes.
   integer                       , intent(inout) :: rc                 !< Return code
 
   ! Local Variables
@@ -88,7 +101,10 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
   real(ESMF_KIND_R8), allocatable :: tauy(:,:)
   real(ESMF_KIND_R8), allocatable :: stkx(:,:,:)
   real(ESMF_KIND_R8), allocatable :: stky(:,:,:)
+  real(ESMF_KIND_R8), allocatable :: work(:,:)
   character(len=*)  , parameter   :: subname = '(mom_import)'
+  character(len=256)              :: stdname
+  integer                         :: field_index
 
   rc = ESMF_SUCCESS
 
@@ -364,6 +380,48 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
       deallocate(stkx,stky)
   endif
 
+  !---
+  ! Tracer flux fields for generic tracers
+  !---
+#ifdef _USE_GENERIC_TRACER
+  if (present(atm_fields)) then
+    ! Set fields in atm_fields from coupler
+    allocate (work(isc:iec,jsc:jec))
+    do n = 1, atm_fields%num_bcs
+      if (atm_fields%bc(n)%flux_type .eq. 'air_sea_deposition') then
+        field_index = ind_deposition
+      elseif (atm_fields%bc(n)%flux_type .eq. 'land_sea_runoff') then
+        field_index = ind_runoff
+      else
+        ! This is a gas flux - set ind_u10 and ind_psurf
+        ! Note, we set these fields even though the pcair field may not be set below. This
+        ! is to allow flux calculation with overridden pcair fields
+        field_index = ind_pcair
+        call set_coupler_type_data(sqrt(ice_ocean_boundary%u10_sqr), n, atm_fields, &
+              idim=(/isc,isc,iec,iec/), jdim=(/jsc,jsc,jec,jec/), field_index=ind_u10)
+        call set_coupler_type_data(ice_ocean_boundary%p, n, atm_fields, &
+              idim=(/isc,isc,iec,iec/), jdim=(/jsc,jsc,jec,jec/), field_index=ind_psurf)
+      endif
+
+      stdname = get_coupled_field_name(atm_fields%bc(n)%name)
+      if (stdname /= UNKNOWN_CMEPS_FIELD) then
+        call ESMF_LogWrite(trim(subname)//': generic_tracer flux, '//trim(atm_fields%bc(n)%name)//&
+          ': setting field index '//CHAR(48+field_index)//' to '//trim(stdname)//' if provided '//&
+          'by coupler, otherwise defaulting to zero', ESMF_LOGMSG_INFO)
+        work(:,:) = 0._ESMF_KIND_R8
+        call state_getimport(importState, trim(stdname), isc, iec, jsc, jec, work, rc=rc)
+        if (ChkErr(rc,__LINE__,u_FILE_u)) return
+        call set_coupler_type_data(work, n, atm_fields, &
+              idim=(/isc,isc,iec,iec/), jdim=(/jsc,jsc,jec,jec/), field_index=field_index)
+      else
+        call ESMF_LogWrite(trim(subname)//': generic_tracer flux, '//trim(atm_fields%bc(n)%name)//&
+              ': no fields set from coupler', ESMF_LOGMSG_INFO)
+      endif
+    enddo
+    deallocate(work)
+  endif
+#endif
+
 end subroutine mom_import
 
 !> Maps outgoing ocean data to ESMF State
diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
index 9ac40daaa4..b6cc4de0be 100644
--- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
+++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
@@ -20,6 +20,7 @@ module MOM_ocean_model_nuopc
 use MOM_constants,           only : CELSIUS_KELVIN_OFFSET, hlf
 use MOM_diag_mediator,       only : diag_ctrl, enable_averages, disable_averaging
 use MOM_diag_mediator,       only : diag_mediator_close_registration, diag_mediator_end
+use MOM_domains,             only : MOM_domain_type, domain2d, clone_MOM_domain, get_domain_extent
 use MOM_domains,             only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE
 use MOM_domains,             only : TO_ALL, Omit_Corners
 use MOM_error_handler,       only : MOM_error, FATAL, WARNING, is_root_pe
@@ -51,8 +52,6 @@ module MOM_ocean_model_nuopc
 use MOM_coupler_types,       only : coupler_type_spawn, coupler_type_write_chksums
 use MOM_coupler_types,       only : coupler_type_initialized, coupler_type_copy_data
 use MOM_coupler_types,       only : coupler_type_set_diags, coupler_type_send_data
-use mpp_domains_mod,         only : domain2d, mpp_get_layout, mpp_get_global_domain
-use mpp_domains_mod,         only : mpp_define_domains, mpp_get_compute_domain, mpp_get_data_domain
 use MOM_io,                  only : stdout
 use MOM_EOS,                 only : gsw_sp_from_sr, gsw_pt_from_ct
 use MOM_wave_interface,      only : wave_parameters_CS, MOM_wave_interface_init
@@ -102,7 +101,7 @@ module MOM_ocean_model_nuopc
                     !! points of the two velocity components. Valid entries
                     !! include AGRID, BGRID_NE, CGRID_NE, BGRID_SW, and CGRID_SW,
                     !! corresponding to the community-standard Arakawa notation.
-                    !! (These are named integers taken from mpp_parameter_mod.)
+                    !! (These are named integers taken from the MOM_domains module.)
                     !! Following MOM5, stagger is BGRID_NE by default when the
                     !! ocean is initialized, but here it is set to -999 so that
                     !! a global max across ocean and non-ocean processors can be
@@ -407,14 +406,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
   ! it also initializes statistical waves.
   call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag)
 
-  if (associated(OS%grid%Domain%maskmap)) then
-    call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, &
-                                      OS%diag, maskmap=OS%grid%Domain%maskmap, &
-                                      gas_fields_ocn=gas_fields_ocn)
-  else
-    call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, &
-                                      OS%diag, gas_fields_ocn=gas_fields_ocn)
-  endif
+  call initialize_ocean_public_type(OS%grid%Domain, Ocean_sfc, OS%diag, gas_fields_ocn=gas_fields_ocn)
 
   ! This call can only occur here if the coupler_bc_type variables have been
   ! initialized already using the information from gas_fields_ocn.
@@ -529,8 +521,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
                           (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.)
 
   ! Translate Ice_ocean_boundary into fluxes.
-  call mpp_get_compute_domain(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), &
-                              index_bnds(3), index_bnds(4))
+  call get_domain_extent(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), index_bnds(3), index_bnds(4))
 
   weight = 1.0
 
@@ -792,7 +783,7 @@ end subroutine ocean_model_end
 subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix)
   type(ocean_state_type),     pointer    :: OS  !< A pointer to the structure containing the
                                                 !! internal ocean state (in).
-  type(time_type),            intent(in) :: Time !< The model time at this call, needed for mpp_write calls.
+  type(time_type),            intent(in) :: Time !< The model time at this call, needed for writing files.
   character(len=*), optional, intent(in) :: directory  !<  An optional directory into which to
                                                 !! write these restart files.
   character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a time-stamp)
@@ -823,16 +814,12 @@ subroutine ocean_model_save_restart(OS, Time, directory, filename_suffix)
 end subroutine ocean_model_save_restart
 
 !> Initialize the public ocean type
-subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, &
-                                        gas_fields_ocn)
-  type(domain2D),          intent(in)    :: input_domain !< The ocean model domain description
+subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, gas_fields_ocn)
+  type(MOM_domain_type),   intent(in)    :: input_domain !< The ocean model domain description
   type(ocean_public_type), intent(inout) :: Ocean_sfc !< A structure containing various publicly
-                                !! visible ocean surface properties after initialization, whose
-                                !! elements are allocated here.
-  type(diag_ctrl),         intent(in)    :: diag  !< A structure that regulates diagnsotic output
-  logical, dimension(:,:), &
-                 optional, intent(in)    :: maskmap !< A mask indicating which virtual processors
-                                              !! are actually in use.  If missing, all are used.
+                                              !! visible ocean surface properties after
+                                              !! initialization, whose elements are allocated here.
+  type(diag_ctrl),         intent(in)    :: diag  !< A structure that regulates diagnostic output
   type(coupler_1d_bc_type), &
                  optional, intent(in)    :: gas_fields_ocn !< If present, this type describes the
                                               !! ocean and surface-ice fields that will participate
@@ -844,14 +831,9 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap,
   ! and have no halos.
   integer :: isc, iec, jsc, jec
 
-  call mpp_get_layout(input_domain,layout)
-  call mpp_get_global_domain(input_domain, xsize=xsz, ysize=ysz)
-  if (PRESENT(maskmap)) then
-    call mpp_define_domains((/1,xsz,1,ysz/),layout,Ocean_sfc%Domain, maskmap=maskmap)
-  else
-    call mpp_define_domains((/1,xsz,1,ysz/),layout,Ocean_sfc%Domain)
-  endif
-  call mpp_get_compute_domain(Ocean_sfc%Domain, isc, iec, jsc, jec)
+  call clone_MOM_domain(input_domain, Ocean_sfc%Domain, halo_size=0, symmetric=.false.)
+
+  call get_domain_extent(Ocean_sfc%Domain, isc, iec, jsc, jec)
 
   allocate ( Ocean_sfc%t_surf (isc:iec,jsc:jec), &
              Ocean_sfc%s_surf (isc:iec,jsc:jec), &
@@ -907,8 +889,7 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   call pass_vector(sfc_state%u, sfc_state%v, G%Domain)
 
-  call mpp_get_compute_domain(Ocean_sfc%Domain, isc_bnd, iec_bnd, &
-                              jsc_bnd, jec_bnd)
+  call get_domain_extent(Ocean_sfc%Domain, isc_bnd, iec_bnd, jsc_bnd, jec_bnd)
   if (present(patm)) then
     ! Check that the inidicies in patm are (isc_bnd:iec_bnd,jsc_bnd:jec_bnd).
     if (.not.present(press_to_z)) call MOM_error(FATAL, &
diff --git a/config_src/external/GFDL_ocean_BGC/generic_tracer.F90 b/config_src/external/GFDL_ocean_BGC/generic_tracer.F90
index 42c386497a..90e7d795ff 100644
--- a/config_src/external/GFDL_ocean_BGC/generic_tracer.F90
+++ b/config_src/external/GFDL_ocean_BGC/generic_tracer.F90
@@ -21,6 +21,7 @@ module generic_tracer
   public generic_tracer_vertdiff_G
   public generic_tracer_get_diag_list
   public generic_tracer_coupler_accumulate
+  public generic_tracer_update_from_coupler
 
   !> Turn on generic tracers (note dangerous use of module data)
   logical :: do_generic_tracer = .true.
@@ -65,6 +66,14 @@ subroutine  generic_tracer_coupler_accumulate(IOB_struc, weight, model_time)
     type(time_type), optional,intent(in) :: model_time !< Time
   end subroutine generic_tracer_coupler_accumulate
 
+  !> Modify the values obtained from the coupler
+  subroutine generic_tracer_update_from_coupler(ilb, jlb, salt_flux_added)
+    integer,                        intent(in) :: ilb    !< Lower bounds of x extent of input arrays on data domain
+    integer,                        intent(in) :: jlb    !< Lower bounds of y extent of input arrays on data domain
+    real, dimension(ilb:,jlb:),     intent(in) :: salt_flux_added !< Surface salt flux into ocean from restoring
+                                                                  !! or flux adjustment [g/m^2/sec]
+  end subroutine generic_tracer_update_from_coupler
+
   !> Calls the corresponding generic_X_update_from_source routine for each package X
   subroutine generic_tracer_source(Temp,Salt,rho_dzt,dzt,hblt_depth,ilb,jlb,tau,dtts,&
        grid_dat,model_time,nbands,max_wavelength_band,sw_pen_band,opacity_band,internal_heat,&
diff --git a/config_src/infra/FMS1/MOM_couplertype_infra.F90 b/config_src/infra/FMS1/MOM_couplertype_infra.F90
index 637f2b5ebf..0fb57c991a 100644
--- a/config_src/infra/FMS1/MOM_couplertype_infra.F90
+++ b/config_src/infra/FMS1/MOM_couplertype_infra.F90
@@ -9,7 +9,10 @@ module MOM_couplertype_infra
 use coupler_types_mod, only : coupler_type_increment_data, coupler_type_rescale_data
 use coupler_types_mod, only : coupler_type_extract_data, coupler_type_set_data
 use coupler_types_mod, only : coupler_type_data_override
-use coupler_types_mod, only : ind_flux, ind_alpha, ind_csurf
+use coupler_types_mod, only : ind_flux, ind_deltap, ind_kw, ind_flux0
+use coupler_types_mod, only : ind_pcair, ind_u10, ind_psurf
+use coupler_types_mod, only : ind_alpha, ind_csurf, ind_sc_no
+use coupler_types_mod, only : ind_runoff, ind_deposition
 use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type
 use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux
 use MOM_domain_infra,  only : domain2D
@@ -22,7 +25,10 @@ module MOM_couplertype_infra
 public :: CT_set_data,  CT_increment_data, CT_rescale_data
 public :: CT_copy_data, CT_extract_data, CT_redistribute_data
 public :: atmos_ocn_coupler_flux
-public :: ind_flux, ind_alpha, ind_csurf
+public :: ind_flux, ind_deltap, ind_kw, ind_flux0
+public :: ind_pcair, ind_u10, ind_psurf
+public :: ind_alpha, ind_csurf, ind_sc_no
+public :: ind_runoff, ind_deposition
 public :: coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type
 
 !> This is the interface to spawn one coupler_bc_type into another.
diff --git a/config_src/infra/FMS2/MOM_couplertype_infra.F90 b/config_src/infra/FMS2/MOM_couplertype_infra.F90
index 3bcccc1dc7..50fd5b09cd 100644
--- a/config_src/infra/FMS2/MOM_couplertype_infra.F90
+++ b/config_src/infra/FMS2/MOM_couplertype_infra.F90
@@ -9,7 +9,10 @@ module MOM_couplertype_infra
 use coupler_types_mod, only : coupler_type_increment_data, coupler_type_rescale_data
 use coupler_types_mod, only : coupler_type_extract_data, coupler_type_set_data
 use coupler_types_mod, only : coupler_type_data_override
-use coupler_types_mod, only : ind_flux, ind_alpha, ind_csurf
+use coupler_types_mod, only : ind_flux, ind_deltap, ind_kw, ind_flux0
+use coupler_types_mod, only : ind_pcair, ind_u10, ind_psurf
+use coupler_types_mod, only : ind_alpha, ind_csurf, ind_sc_no
+use coupler_types_mod, only : ind_runoff, ind_deposition
 use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type
 use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux
 use MOM_domain_infra,  only : domain2D
@@ -22,7 +25,10 @@ module MOM_couplertype_infra
 public :: CT_set_data,  CT_increment_data, CT_rescale_data
 public :: CT_copy_data, CT_extract_data, CT_redistribute_data
 public :: atmos_ocn_coupler_flux
-public :: ind_flux, ind_alpha, ind_csurf
+public :: ind_flux, ind_deltap, ind_kw, ind_flux0
+public :: ind_pcair, ind_u10, ind_psurf
+public :: ind_alpha, ind_csurf, ind_sc_no
+public :: ind_runoff, ind_deposition
 public :: coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type
 
 !> This is the interface to spawn one coupler_bc_type into another.
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index 998713d1f1..a2efcaac7d 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -2280,6 +2280,11 @@ subroutine fluxes_accumulate(flux_tmp, fluxes, G, wt2, forces)
 
     fluxes%salt_flux(i,j) = wt1*fluxes%salt_flux(i,j) + wt2*flux_tmp%salt_flux(i,j)
   enddo ; enddo
+  if (associated(fluxes%salt_flux_added) .and. associated(flux_tmp%salt_flux_added)) then
+    do j=js,je ; do i=is,ie
+      fluxes%salt_flux_added(i,j) = wt1*fluxes%salt_flux_added(i,j) + wt2*flux_tmp%salt_flux_added(i,j)
+    enddo ; enddo
+  endif
   if (associated(fluxes%heat_added) .and. associated(flux_tmp%heat_added)) then
     do j=js,je ; do i=is,ie
       fluxes%heat_added(i,j) = wt1*fluxes%heat_added(i,j) + wt2*flux_tmp%heat_added(i,j)
diff --git a/src/framework/MOM_coupler_types.F90 b/src/framework/MOM_coupler_types.F90
index b931a2ddd2..cac9309b5e 100644
--- a/src/framework/MOM_coupler_types.F90
+++ b/src/framework/MOM_coupler_types.F90
@@ -9,7 +9,10 @@ module MOM_coupler_types
 use MOM_couplertype_infra, only : CT_copy_data, CT_increment_data, CT_rescale_data
 use MOM_couplertype_infra, only : CT_set_data, CT_extract_data, CT_redistribute_data
 use MOM_couplertype_infra, only : coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type
-use MOM_couplertype_infra, only : ind_flux, ind_alpha, ind_csurf
+use MOM_couplertype_infra, only : ind_flux, ind_deltap, ind_kw, ind_flux0
+use MOM_couplertype_infra, only : ind_pcair, ind_u10, ind_psurf
+use MOM_couplertype_infra, only : ind_alpha, ind_csurf, ind_sc_no
+use MOM_couplertype_infra, only : ind_runoff, ind_deposition
 use MOM_domain_infra,      only : domain2D
 use MOM_time_manager,      only : time_type
 
@@ -23,7 +26,10 @@ module MOM_coupler_types
 public :: coupler_1d_bc_type, coupler_2d_bc_type, coupler_3d_bc_type
 ! These are encoding constant parameters that indicate whether a flux, solubility or
 ! surface ocean concentration are being set or accessed with an inquiry.
-public :: ind_flux, ind_alpha, ind_csurf
+public :: ind_flux, ind_deltap, ind_kw, ind_flux0
+public :: ind_pcair, ind_u10, ind_psurf
+public :: ind_alpha, ind_csurf, ind_sc_no
+public :: ind_runoff, ind_deposition
 
 !> This is the interface to spawn one coupler_bc_type into another.
 interface coupler_type_spawn
diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90
index 6b10d15e2f..e092aa7d6d 100644
--- a/src/tracer/MOM_generic_tracer.F90
+++ b/src/tracer/MOM_generic_tracer.F90
@@ -21,7 +21,7 @@ module MOM_generic_tracer
   use generic_tracer, only: generic_tracer_coupler_get, generic_tracer_coupler_set
   use generic_tracer, only: generic_tracer_end, generic_tracer_get_list, do_generic_tracer
   use generic_tracer, only: generic_tracer_update_from_bottom,generic_tracer_vertdiff_G
-  use generic_tracer, only: generic_tracer_coupler_accumulate
+  use generic_tracer, only: generic_tracer_coupler_accumulate, generic_tracer_update_from_coupler
 
   use g_tracer_utils,   only: g_tracer_get_name,g_tracer_set_values,g_tracer_set_common,g_tracer_get_common
   use g_tracer_utils,   only: g_tracer_get_next,g_tracer_type,g_tracer_is_prog,g_tracer_flux_init
@@ -527,6 +527,11 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
     !      the fluxes without coming into this subroutine.
     !      MOM5 has to modified to conform.
 
+    !
+    !Call the generic_tracer's update_from_coupler routine (convert salt_flux_added to g/m^2/sec)
+    !
+    call generic_tracer_update_from_coupler(G%isd, G%jsd, 1000*(US%RZ_T_to_kg_m2s*fluxes%salt_flux_added))
+
     !
     !Add contribution of river to surface flux
     !