MOM_OCMIP2_CFC.F90

!> Simulates CFCs using the OCMIP2 protocols
module mom_ocmip2_cfc

! This file is part of MOM6. See LICENSE.md for the license.

use mom_diag_mediator, only : diag_ctrl
use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_hor_index, only : hor_index_type
use mom_grid, only : ocean_grid_type
use mom_io, only : file_exists, mom_read_data, slasher, vardesc, var_desc, query_vardesc
use mom_open_boundary, only : ocean_obc_type
use mom_restart, only : query_initialized, mom_restart_cs
use mom_sponge, only : set_up_sponge_field, sponge_cs
use mom_time_manager, only : time_type
use mom_tracer_registry, only : register_tracer, tracer_registry_type
use mom_tracer_diabatic, only : tracer_vertdiff, applytracerboundaryfluxesinout
use mom_tracer_z_init, only : tracer_z_init
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : surface
use mom_verticalgrid, only : verticalgrid_type

use coupler_types_mod, only : ind_flux, ind_alpha, ind_csurf
use coupler_types_mod, only : coupler_type_extract_data, coupler_type_set_data
use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux

implicit none ; private

#include <MOM_memory.h>

public register_ocmip2_cfc, initialize_ocmip2_cfc, flux_init_ocmip2_cfc
public ocmip2_cfc_column_physics, ocmip2_cfc_surface_state
public ocmip2_cfc_stock, ocmip2_cfc_end


integer, parameter :: ntr = 2 !< the number of tracers in this module.

!> The control structure for the  OCMPI2_CFC tracer package
type, public :: ocmip2_cfc_cs ; private
  character(len=200) :: ic_file !< The file in which the CFC initial values can
                                !! be found, or an empty string for internal initilaization.
  logical :: z_ic_file !< If true, the IC_file is in Z-space.  The default is false..
  type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
  type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the MOM6 tracer registry
  real, pointer, dimension(:,:,:) :: &
    cfc11 => null(), &     !< The CFC11 concentration [mol m-3].
    cfc12 => null()        !< The CFC12 concentration [mol m-3].
  ! In the following variables a suffix of _11 refers to CFC11 and _12 to CFC12.
  !>@{ Coefficients used in the CFC11 and CFC12 solubility calculation
  real :: a1_11, a1_12   ! Coefficients for calculating CFC11 and CFC12 Schmidt numbers [nondim]
  real :: a2_11, a2_12   ! Coefficients for calculating CFC11 and CFC12 Schmidt numbers [degC-1]
  real :: a3_11, a3_12   ! Coefficients for calculating CFC11 and CFC12 Schmidt numbers [degC-2]
  real :: a4_11, a4_12   ! Coefficients for calculating CFC11 and CFC12 Schmidt numbers [degC-3]

  real :: d1_11, d1_12   ! Coefficients for calculating CFC11 and CFC12 solubilities [nondim]
  real :: d2_11, d2_12   ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-1]
  real :: d3_11, d3_12   ! Coefficients for calculating CFC11 and CFC12 solubilities [log(hectoKelvin)-1]
  real :: d4_11, d4_12   ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-2]

  real :: e1_11, e1_12   ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1]
  real :: e2_11, e2_12   ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1 hectoKelvin-1]
  real :: e3_11, e3_12   ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-2 hectoKelvin-2]
  !>@}
  real :: cfc11_ic_val = 0.0    !< The initial value assigned to CFC11 [mol m-3].
  real :: cfc12_ic_val = 0.0    !< The initial value assigned to CFC12 [mol m-3].
  real :: cfc11_land_val = -1.0 !< The value of CFC11 used where land is masked out [mol m-3].
  real :: cfc12_land_val = -1.0 !< The value of CFC12 used where land is masked out [mol m-3].
  logical :: tracers_may_reinit !< If true, tracers may be reset via the initialization code
                                !! if they are not found in the restart files.
  character(len=16) :: cfc11_name !< CFC11 variable name
  character(len=16) :: cfc12_name !< CFC12 variable name

  integer :: ind_cfc_11_flux  !< Index returned by aof_set_coupler_flux that is used to
                              !! pack and unpack surface boundary condition arrays.
  integer :: ind_cfc_12_flux  !< Index returned by aof_set_coupler_flux that is used to
                              !! pack and unpack surface boundary condition arrays.

  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate
                                             !! the timing of diagnostic output.
  type(mom_restart_cs), pointer :: restart_csp => null()  !< Model restart control structure

  ! The following vardesc types contain a package of metadata about each tracer.
  type(vardesc) :: cfc11_desc !< A set of metadata for the CFC11 tracer
  type(vardesc) :: cfc12_desc !< A set of metadata for the CFC12 tracer
end type ocmip2_cfc_cs

contains

!> Register the OCMIP2 CFC tracers to be used with MOM and read the parameters
!! that are used with this tracer package
function register_ocmip2_cfc(HI, GV, param_file, CS, tr_Reg, restart_CS)
  type(hor_index_type),    intent(in) :: HI         !< A horizontal index type structure.
  type(verticalgrid_type), intent(in) :: GV         !< The ocean's vertical grid structure.
  type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time parameters.
  type(ocmip2_cfc_cs),     pointer    :: CS         !< A pointer that is set to point to the control
                                                    !! structure for this module.
  type(tracer_registry_type), &
                           pointer    :: tr_Reg     !< A pointer to the tracer registry.
  type(mom_restart_cs),    pointer    :: restart_CS !< A pointer to the restart control structure.
! This subroutine is used to register tracer fields and subroutines
! to be used with MOM.

  ! Local variables
  character(len=40)  :: mdl = "MOM_OCMIP2_CFC" ! This module's name.
  character(len=200) :: inputdir ! The directory where NetCDF input files are.
  ! This include declares and sets the variable "version".
#include "version_variable.h"
  real, dimension(:,:,:), pointer :: tr_ptr => null()
  real :: a11_dflt(4), a12_dflt(4) ! Default values of the various coefficients
  real :: d11_dflt(4), d12_dflt(4) ! In the expressions for the solubility and
  real :: e11_dflt(3), e12_dflt(3) ! Schmidt numbers.
  character(len=48) :: flux_units ! The units for tracer fluxes.
  logical :: register_OCMIP2_CFC
  integer :: isd, ied, jsd, jed, nz, m

  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke

  if (associated(cs)) then
    call mom_error(warning, "register_OCMIP2_CFC called with an "// &
                            "associated control structure.")
    return
  endif
  allocate(cs)

  ! This call sets default properties for the air-sea CFC fluxes and obtains the
  ! indicies for the CFC11 and CFC12 flux coupling.
  call flux_init_ocmip2_cfc(cs, verbosity=3)
  if ((cs%ind_cfc_11_flux < 0) .or. (cs%ind_cfc_11_flux < 0)) then
    ! This is most likely to happen with the dummy version of aof_set_coupler_flux
    ! used in ocean-only runs.
    call mom_error(warning, "CFCs are currently only set up to be run in " // &
                   " coupled model configurations, and will be disabled.")
    deallocate(cs)
    register_ocmip2_cfc = .false.
    return
  endif

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "CFC_IC_FILE", cs%IC_file, &
                 "The file in which the CFC initial values can be "//&
                 "found, or an empty string for internal initialization.", &
                 default=" ")
  if ((len_trim(cs%IC_file) > 0) .and. (scan(cs%IC_file,'/') == 0)) then
    ! Add the directory if CS%IC_file is not already a complete path.
    call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
    cs%IC_file = trim(slasher(inputdir))//trim(cs%IC_file)
    call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", cs%IC_file)
  endif
  call get_param(param_file, mdl, "CFC_IC_FILE_IS_Z", cs%Z_IC_file, &
                 "If true, CFC_IC_FILE is in depth space, not layer space", &
                 default=.false.)
  call get_param(param_file, mdl, "TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
                 "If true, tracers may go through the initialization code "//&
                 "if they are not found in the restart files.  Otherwise "//&
                 "it is a fatal error if tracers are not found in the "//&
                 "restart files of a restarted run.", default=.false.)

  !   The following vardesc types contain a package of metadata about each tracer,
  ! including, the name; units; longname; and grid information.
  cs%CFC11_name = "CFC11" ; cs%CFC12_name = "CFC12"
  cs%CFC11_desc = var_desc(cs%CFC11_name,"mol m-3","CFC-11 Concentration", caller=mdl)
  cs%CFC12_desc = var_desc(cs%CFC12_name,"mol m-3","CFC-12 Concentration", caller=mdl)
  if (gv%Boussinesq) then ; flux_units = "mol s-1"
  else ; flux_units = "mol m-3 kg s-1" ; endif

  allocate(cs%CFC11(isd:ied,jsd:jed,nz)) ; cs%CFC11(:,:,:) = 0.0
  allocate(cs%CFC12(isd:ied,jsd:jed,nz)) ; cs%CFC12(:,:,:) = 0.0

  ! This pointer assignment is needed to force the compiler not to do a copy in
  ! the registration calls.  Curses on the designers and implementers of F90.
  tr_ptr => cs%CFC11
  ! Register CFC11 for horizontal advection, diffusion, and restarts.
  call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
                       tr_desc=cs%CFC11_desc, registry_diags=.true., &
                       flux_units=flux_units, &
                       restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)
  ! Do the same for CFC12
  tr_ptr => cs%CFC12
  call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
                       tr_desc=cs%CFC12_desc, registry_diags=.true., &
                       flux_units=flux_units, &
                       restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)

  ! Set and read the various empirical coefficients.

!-----------------------------------------------------------------------
! Default Schmidt number coefficients for CFC11 (_11) and CFC12 (_12) are given
! by Zheng et al (1998), JGR vol 103, C1.
!-----------------------------------------------------------------------
  a11_dflt(:) = (/ 3501.8, -210.31,  6.1851, -0.07513 /)
  a12_dflt(:) = (/ 3845.4, -228.95,  6.1908, -0.06743 /)
  call get_param(param_file, mdl, "CFC11_A1", cs%a1_11, &
                 "A coefficient in the Schmidt number of CFC11.", &
                 units="nondim", default=a11_dflt(1))
  call get_param(param_file, mdl, "CFC11_A2", cs%a2_11, &
                 "A coefficient in the Schmidt number of CFC11.", &
                 units="degC-1", default=a11_dflt(2))
  call get_param(param_file, mdl, "CFC11_A3", cs%a3_11, &
                 "A coefficient in the Schmidt number of CFC11.", &
                 units="degC-2", default=a11_dflt(3))
  call get_param(param_file, mdl, "CFC11_A4", cs%a4_11, &
                 "A coefficient in the Schmidt number of CFC11.", &
                 units="degC-3", default=a11_dflt(4))

  call get_param(param_file, mdl, "CFC12_A1", cs%a1_12, &
                 "A coefficient in the Schmidt number of CFC12.", &
                 units="nondim", default=a12_dflt(1))
  call get_param(param_file, mdl, "CFC12_A2", cs%a2_12, &
                 "A coefficient in the Schmidt number of CFC12.", &
                 units="degC-1", default=a12_dflt(2))
  call get_param(param_file, mdl, "CFC12_A3", cs%a3_12, &
                 "A coefficient in the Schmidt number of CFC12.", &
                 units="degC-2", default=a12_dflt(3))
  call get_param(param_file, mdl, "CFC12_A4", cs%a4_12, &
                 "A coefficient in the Schmidt number of CFC12.", &
                 units="degC-3", default=a12_dflt(4))

!-----------------------------------------------------------------------
! Solubility coefficients for alpha in mol/l/atm for CFC11 (_11) and CFC12 (_12)
! after Warner and Weiss (1985) DSR, vol 32.
!-----------------------------------------------------------------------
  d11_dflt(:) = (/ -229.9261, 319.6552, 119.4471, -1.39165 /)
  e11_dflt(:) = (/ -0.142382, 0.091459, -0.0157274 /)
  d12_dflt(:) = (/ -218.0971, 298.9702, 113.8049, -1.39165 /)
  e12_dflt(:) = (/ -0.143566, 0.091015, -0.0153924 /)

  call get_param(param_file, mdl, "CFC11_D1", cs%d1_11, &
                 "A coefficient in the solubility of CFC11.", &
                 units="none", default=d11_dflt(1))
  call get_param(param_file, mdl, "CFC11_D2", cs%d2_11, &
                 "A coefficient in the solubility of CFC11.", &
                 units="hK", default=d11_dflt(2))
  call get_param(param_file, mdl, "CFC11_D3", cs%d3_11, &
                 "A coefficient in the solubility of CFC11.", &
                 units="none", default=d11_dflt(3))
  call get_param(param_file, mdl, "CFC11_D4", cs%d4_11, &
                 "A coefficient in the solubility of CFC11.", &
                 units="hK-2", default=d11_dflt(4))
  call get_param(param_file, mdl, "CFC11_E1", cs%e1_11, &
                 "A coefficient in the solubility of CFC11.", &
                 units="PSU-1", default=e11_dflt(1))
  call get_param(param_file, mdl, "CFC11_E2", cs%e2_11, &
                 "A coefficient in the solubility of CFC11.", &
                 units="PSU-1 hK-1", default=e11_dflt(2))
  call get_param(param_file, mdl, "CFC11_E3", cs%e3_11, &
                 "A coefficient in the solubility of CFC11.", &
                 units="PSU-1 hK-2", default=e11_dflt(3))

  call get_param(param_file, mdl, "CFC12_D1", cs%d1_12, &
                 "A coefficient in the solubility of CFC12.", &
                 units="none", default=d12_dflt(1))
  call get_param(param_file, mdl, "CFC12_D2", cs%d2_12, &
                 "A coefficient in the solubility of CFC12.", &
                 units="hK", default=d12_dflt(2))
  call get_param(param_file, mdl, "CFC12_D3", cs%d3_12, &
                 "A coefficient in the solubility of CFC12.", &
                 units="none", default=d12_dflt(3))
  call get_param(param_file, mdl, "CFC12_D4", cs%d4_12, &
                 "A coefficient in the solubility of CFC12.", &
                 units="hK-2", default=d12_dflt(4))
  call get_param(param_file, mdl, "CFC12_E1", cs%e1_12, &
                 "A coefficient in the solubility of CFC12.", &
                 units="PSU-1", default=e12_dflt(1))
  call get_param(param_file, mdl, "CFC12_E2", cs%e2_12, &
                 "A coefficient in the solubility of CFC12.", &
                 units="PSU-1 hK-1", default=e12_dflt(2))
  call get_param(param_file, mdl, "CFC12_E3", cs%e3_12, &
                 "A coefficient in the solubility of CFC12.", &
                 units="PSU-1 hK-2", default=e12_dflt(3))

  cs%tr_Reg => tr_reg
  cs%restart_CSp => restart_cs

  register_ocmip2_cfc = .true.
end function register_ocmip2_cfc

!> This subroutine initializes the air-sea CFC fluxes, and optionally returns
!! the indicies of these fluxes.  It can safely be called multiple times.
subroutine flux_init_ocmip2_cfc(CS, verbosity)
  type(ocmip2_cfc_cs), optional, pointer :: CS !< An optional pointer to the control structure
                                               !! for this module; if not present, the flux indicies
                                               !! are not stored.
  integer,             optional, intent(in) :: verbosity !< A 0-9 integer indicating a level of verbosity.

  ! These can be overridden later in via the field manager?
  character(len=128) :: default_ice_restart_file = 'ice_ocmip2_cfc.res.nc'
  character(len=128) :: default_ocean_restart_file = 'ocmip2_cfc.res.nc'
  integer :: ind_flux(2) ! Integer indices of the fluxes

  ! These calls obtain the indices for the CFC11 and CFC12 flux coupling.  They
  ! can safely be called multiple times.
  ind_flux(1) = aof_set_coupler_flux('cfc_11_flux', &
       flux_type = 'air_sea_gas_flux', implementation = 'ocmip2', &
       param = (/ 9.36e-07, 9.7561e-06 /), &
       ice_restart_file = default_ice_restart_file, &
       ocean_restart_file = default_ocean_restart_file, &
       caller = "register_OCMIP2_CFC", verbosity=verbosity)
  ind_flux(2) = aof_set_coupler_flux('cfc_12_flux', &
       flux_type = 'air_sea_gas_flux', implementation = 'ocmip2', &
       param = (/ 9.36e-07, 9.7561e-06 /), &
       ice_restart_file = default_ice_restart_file, &
       ocean_restart_file = default_ocean_restart_file, &
       caller = "register_OCMIP2_CFC", verbosity=verbosity)

  if (present(cs)) then ; if (associated(cs)) then
    cs%ind_cfc_11_flux = ind_flux(1)
    cs%ind_cfc_12_flux = ind_flux(2)
  endif ; endif

end subroutine flux_init_ocmip2_cfc

!> Initialize the OCMP2 CFC tracer fields and set up the tracer output.
subroutine initialize_ocmip2_cfc(restart, day, G, GV, US, h, diag, OBC, CS, &
                                 sponge_CSp)
  logical,                        intent(in) :: restart    !< .true. if the fields have already been
                                                           !! read from a restart file.
  type(time_type), target,        intent(in) :: day        !< Time of the start of the run.
  type(ocean_grid_type),          intent(in) :: G          !< The ocean's grid structure.
  type(verticalgrid_type),        intent(in) :: GV         !< The ocean's vertical grid structure.
  type(unit_scale_type),          intent(in) :: US         !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                  intent(in) :: h          !< Layer thicknesses [H ~> m or kg m-2].
  type(diag_ctrl), target,        intent(in) :: diag       !< A structure that is used to regulate
                                                           !! diagnostic output.
  type(ocean_obc_type),           pointer    :: OBC        !< This open boundary condition type
                                                           !! specifies whether, where, and what
                                                           !! open boundary conditions are used.
  type(ocmip2_cfc_cs),            pointer    :: CS         !< The control structure returned by a
                                                           !! previous call to register_OCMIP2_CFC.
  type(sponge_cs),                pointer    :: sponge_CSp !< A pointer to the control structure for
                                                           !! the sponges, if they are in use.
                                                           !! Otherwise this may be unassociated.
!   This subroutine initializes the NTR tracer fields in tr(:,:,:,:)
! and it sets up the tracer output.

  logical :: from_file = .false.

  if (.not.associated(cs)) return

  cs%Time => day
  cs%diag => diag

  if (.not.restart .or. (cs%tracers_may_reinit .and. &
      .not.query_initialized(cs%CFC11, cs%CFC11_name, cs%restart_CSp))) &
    call init_tracer_cfc(h, cs%CFC11, cs%CFC11_name, cs%CFC11_land_val, &
                         cs%CFC11_IC_val, g, us, cs)

  if (.not.restart .or. (cs%tracers_may_reinit .and. &
      .not.query_initialized(cs%CFC12, cs%CFC12_name, cs%restart_CSp))) &
    call init_tracer_cfc(h, cs%CFC12, cs%CFC12_name, cs%CFC12_land_val, &
                         cs%CFC12_IC_val, g, us, cs)

  if (associated(obc)) then
  ! Steal from updated DOME in the fullness of time.
  endif

end subroutine initialize_ocmip2_cfc

!>This subroutine initializes a tracer array.
subroutine init_tracer_cfc(h, tr, name, land_val, IC_val, G, US, CS)
  type(ocean_grid_type),                    intent(in)  :: G    !< The ocean's grid structure
  type(unit_scale_type),                    intent(in)  :: US   !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)  :: h    !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: tr   !< The tracer concentration array
  character(len=*),                         intent(in)  :: name !< The tracer name
  real,                                     intent(in)  :: land_val !< A value the tracer takes over land
  real,                                     intent(in)  :: IC_val !< The initial condition value for the tracer
  type(ocmip2_cfc_cs),                      pointer     :: CS   !< The control structure returned by a
                                                                !! previous call to register_OCMIP2_CFC.

  ! This subroutine initializes a tracer array.

  logical :: OK
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  if (len_trim(cs%IC_file) > 0) then
    !  Read the tracer concentrations from a netcdf file.
    if (.not.file_exists(cs%IC_file, g%Domain)) &
      call mom_error(fatal, "initialize_OCMIP2_CFC: Unable to open "//cs%IC_file)
    if (cs%Z_IC_file) then
      ok = tracer_z_init(tr, h, cs%IC_file, name, g, us)
      if (.not.ok) then
        ok = tracer_z_init(tr, h, cs%IC_file, trim(name), g, us)
        if (.not.ok) call mom_error(fatal,"initialize_OCMIP2_CFC: "//&
                "Unable to read "//trim(name)//" from "//&
                trim(cs%IC_file)//".")
      endif
    else
      call mom_read_data(cs%IC_file, trim(name), tr, g%Domain)
    endif
  else
    do k=1,nz ; do j=js,je ; do i=is,ie
      if (g%mask2dT(i,j) < 0.5) then
        tr(i,j,k) = land_val
      else
        tr(i,j,k) = ic_val
      endif
    enddo ; enddo ; enddo
  endif

end subroutine init_tracer_cfc

!>  This subroutine applies diapycnal diffusion, souces and sinks and any other column
!! tracer physics or chemistry to the OCMIP2 CFC tracers.
!! CFCs are relatively simple, as they are passive tracers with only a surface flux as a source.
subroutine ocmip2_cfc_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, &
              evap_CFL_limit, minimum_forcing_depth)
  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure
  type(verticalgrid_type), intent(in) :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: ea   !< an array to which the amount of fluid entrained
                                              !! from the layer above during this call will be
                                              !! added [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: eb   !< an array to which the amount of fluid entrained
                                              !! from the layer below during this call will be
                                              !! added [H ~> m or kg m-2].
  type(forcing),           intent(in) :: fluxes !< A structure containing pointers to thermodynamic
                                              !! and tracer forcing fields.  Unused fields have NULL ptrs.
  real,                    intent(in) :: dt   !< The amount of time covered by this call [T ~> s]
  type(unit_scale_type),   intent(in) :: US   !< A dimensional unit scaling type
  type(ocmip2_cfc_cs),     pointer    :: CS   !< The control structure returned by a
                                              !! previous call to register_OCMIP2_CFC.
  real,          optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can
                                              !! be fluxed out of the top layer in a timestep [nondim]
  real,          optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
                                              !! fluxes can be applied [H ~> m or kg m-2]
!   This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! CFCs are relatively simple, as they are passive tracers. with only a surface
! flux as a source.

! The arguments to this subroutine are redundant in that
!     h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)

  ! Local variables
  real :: b1(szi_(g))          ! b1 and c1 are variables used by the
  real :: c1(szi_(g),szk_(g))  ! tridiagonal solver.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    CFC11_flux, &    ! The fluxes of CFC11 and CFC12 into the ocean, in the
    CFC12_flux       ! units of CFC concentrations times meters per second.
  real, pointer, dimension(:,:,:) :: CFC11 => null(), cfc12 => null()
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
  integer :: i, j, k, m, is, ie, js, je, nz, idim(4), jdim(4)

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
  idim(:) = (/g%isd, is, ie, g%ied/) ; jdim(:) = (/g%jsd, js, je, g%jed/)

  if (.not.associated(cs)) return

  cfc11 => cs%CFC11 ; cfc12 => cs%CFC12

  ! These two calls unpack the fluxes from the input arrays.
  !   The -GV%Rho0 changes the sign convention of the flux and changes the units
  ! of the flux from [Conc. m s-1] to [Conc. kg m-2 T-1].
  call coupler_type_extract_data(fluxes%tr_fluxes, cs%ind_cfc_11_flux, ind_flux, cfc11_flux, &
                                 scale_factor=-g%US%R_to_kg_m3*gv%Rho0*us%T_to_s, idim=idim, jdim=jdim)
  call coupler_type_extract_data(fluxes%tr_fluxes, cs%ind_cfc_12_flux, ind_flux, cfc12_flux, &
                                 scale_factor=-g%US%R_to_kg_m3*gv%Rho0*us%T_to_s, idim=idim, jdim=jdim)

  ! Use a tridiagonal solver to determine the concentrations after the
  ! surface source is applied and diapycnal advection and diffusion occurs.
  if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
    do k=1,nz ;do j=js,je ; do i=is,ie
      h_work(i,j,k) = h_old(i,j,k)
    enddo ; enddo ; enddo
    call applytracerboundaryfluxesinout(g, gv, cfc11, dt, fluxes, h_work, &
                                        evap_cfl_limit, minimum_forcing_depth)
    call tracer_vertdiff(h_work, ea, eb, dt, cfc11, g, gv, sfc_flux=cfc11_flux)

    do k=1,nz ;do j=js,je ; do i=is,ie
      h_work(i,j,k) = h_old(i,j,k)
    enddo ; enddo ; enddo
    call applytracerboundaryfluxesinout(g, gv, cfc12, dt, fluxes, h_work, &
                                        evap_cfl_limit, minimum_forcing_depth)
    call tracer_vertdiff(h_work, ea, eb, dt, cfc12, g, gv, sfc_flux=cfc12_flux)
  else
    call tracer_vertdiff(h_old, ea, eb, dt, cfc11, g, gv, sfc_flux=cfc11_flux)
    call tracer_vertdiff(h_old, ea, eb, dt, cfc12, g, gv, sfc_flux=cfc12_flux)
  endif

  ! Write out any desired diagnostics from tracer sources & sinks here.

end subroutine ocmip2_cfc_column_physics

!> This function calculates the mass-weighted integral of all tracer stocks,
!! returning the number of stocks it has calculated.  If the stock_index
!! is present, only the stock corresponding to that coded index is returned.
function ocmip2_cfc_stock(h, stocks, G, GV, CS, names, units, stock_index)
  type(ocean_grid_type),           intent(in)    :: G      !< The ocean's grid structure.
  type(verticalgrid_type),         intent(in)    :: GV     !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                   intent(in)    :: h      !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(:),              intent(out)   :: stocks !< the mass-weighted integrated amount of each
                                                           !! tracer, in kg times concentration units [kg conc].
  type(ocmip2_cfc_cs),             pointer       :: CS     !< The control structure returned by a
                                                           !! previous call to register_OCMIP2_CFC.
  character(len=*), dimension(:),  intent(out)   :: names  !< The names of the stocks calculated.
  character(len=*), dimension(:),  intent(out)   :: units  !< The units of the stocks calculated.
  integer, optional,               intent(in)    :: stock_index !< The coded index of a specific
                                                                !! stock being sought.
  integer                                        :: OCMIP2_CFC_stock !< The number of stocks calculated here.

  ! Local variables
  real :: mass
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  ocmip2_cfc_stock = 0
  if (.not.associated(cs)) return

  if (present(stock_index)) then ; if (stock_index > 0) then
    ! Check whether this stock is available from this routine.

    ! No stocks from this routine are being checked yet.  Return 0.
    return
  endif ; endif

  call query_vardesc(cs%CFC11_desc, name=names(1), units=units(1), caller="OCMIP2_CFC_stock")
  call query_vardesc(cs%CFC12_desc, name=names(2), units=units(2), caller="OCMIP2_CFC_stock")
  units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg"

  stocks(1) = 0.0 ; stocks(2) = 0.0
  do k=1,nz ; do j=js,je ; do i=is,ie
    mass = g%mask2dT(i,j) * g%US%L_to_m**2*g%areaT(i,j) * h(i,j,k)
    stocks(1) = stocks(1) + cs%CFC11(i,j,k) * mass
    stocks(2) = stocks(2) + cs%CFC12(i,j,k) * mass
  enddo ; enddo ; enddo
  stocks(1) = gv%H_to_kg_m2 * stocks(1)
  stocks(2) = gv%H_to_kg_m2 * stocks(2)

  ocmip2_cfc_stock = 2

end function ocmip2_cfc_stock

!> This subroutine extracts the surface CFC concentrations and other fields that
!! are shared with the atmosphere to calculate CFC fluxes.
subroutine ocmip2_cfc_surface_state(sfc_state, h, G, CS)
  type(ocean_grid_type),  intent(in)    :: G  !< The ocean's grid structure.
  type(surface),          intent(inout) :: sfc_state !< A structure containing fields that
                                              !! describe the surface state of the ocean.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                          intent(in)    :: h  !< Layer thickness [H ~> m or kg m-2].
  type(ocmip2_cfc_cs),    pointer       :: CS !< The control structure returned by a previous
                                              !! call to register_OCMIP2_CFC.

  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: &
    CFC11_Csurf, &  ! The CFC-11 surface concentrations times the Schmidt number term [mol m-3].
    CFC12_Csurf, &  ! The CFC-12 surface concentrations times the Schmidt number term [mol m-3].
    CFC11_alpha, &  ! The CFC-11 solubility [mol m-3 pptv-1].
    CFC12_alpha     ! The CFC-12 solubility [mol m-3 pptv-1].
  real :: ta        ! Absolute sea surface temperature [hectoKelvin] (Why use such bizzare units?)
  real :: sal       ! Surface salinity [PSU].
  real :: SST       ! Sea surface temperature [degC].
  real :: alpha_11  ! The solubility of CFC 11 [mol m-3 pptv-1].
  real :: alpha_12  ! The solubility of CFC 12 [mol m-3 pptv-1].
  real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12.
  real :: sc_no_term   ! A term related to the Schmidt number.
  integer :: i, j, m, is, ie, js, je, idim(4), jdim(4)

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  idim(:) = (/g%isd, is, ie, g%ied/) ; jdim(:) = (/g%jsd, js, je, g%jed/)

  if (.not.associated(cs)) return

  do j=js,je ; do i=is,ie
    ta = max(0.01, (sfc_state%SST(i,j) + 273.15) * 0.01) ! Why is this in hectoKelvin?
    sal = sfc_state%SSS(i,j) ; sst = sfc_state%SST(i,j)
    !    Calculate solubilities using Warner and Weiss (1985) DSR, vol 32.
    ! The final result is in mol/cm3/pptv (1 part per trillion 1e-12)
    ! Use Bullister and Wisegavger for CCl4.
    ! The factor 1.e-09 converts from mol/(l * atm) to mol/(m3 * pptv).
    alpha_11 = exp(cs%d1_11 + cs%d2_11/ta + cs%d3_11*log(ta) + cs%d4_11*ta**2 +&
                   sal * ((cs%e3_11 * ta + cs%e2_11) * ta + cs%e1_11)) * &
               1.0e-09 * g%mask2dT(i,j)
    alpha_12 = exp(cs%d1_12 + cs%d2_12/ta + cs%d3_12*log(ta) + cs%d4_12*ta**2 +&
                   sal * ((cs%e3_12 * ta + cs%e2_12) * ta + cs%e1_12)) * &
               1.0e-09 * g%mask2dT(i,j)
    !   Calculate Schmidt numbers using coefficients given by
    ! Zheng et al (1998), JGR vol 103, C1.
    sc_11 = cs%a1_11 + sst * (cs%a2_11 + sst * (cs%a3_11 + sst * cs%a4_11)) * &
            g%mask2dT(i,j)
    sc_12 = cs%a1_12 + sst * (cs%a2_12 + sst * (cs%a3_12 + sst * cs%a4_12)) * &
            g%mask2dT(i,j)
    ! The abs here is to avoid NaNs. The model should be failing at this point.
    sc_no_term = sqrt(660.0 / (abs(sc_11) + 1.0e-30))
    cfc11_alpha(i,j) = alpha_11 * sc_no_term
    cfc11_csurf(i,j) = cs%CFC11(i,j,1) * sc_no_term

    sc_no_term = sqrt(660.0 / (abs(sc_12) + 1.0e-30))
    cfc12_alpha(i,j) = alpha_12 * sc_no_term
    cfc12_csurf(i,j) = cs%CFC12(i,j,1) * sc_no_term
  enddo ; enddo

  !   These calls load these values into the appropriate arrays in the
  ! coupler-type structure.
  call coupler_type_set_data(cfc11_alpha, cs%ind_cfc_11_flux, ind_alpha, &
                             sfc_state%tr_fields, idim=idim, jdim=jdim)
  call coupler_type_set_data(cfc11_csurf, cs%ind_cfc_11_flux, ind_csurf, &
                             sfc_state%tr_fields, idim=idim, jdim=jdim)
  call coupler_type_set_data(cfc12_alpha, cs%ind_cfc_12_flux, ind_alpha, &
                             sfc_state%tr_fields, idim=idim, jdim=jdim)
  call coupler_type_set_data(cfc12_csurf, cs%ind_cfc_12_flux, ind_csurf, &
                             sfc_state%tr_fields, idim=idim, jdim=jdim)

end subroutine ocmip2_cfc_surface_state

!> Deallocate any memory associated with the OCMIP2 CFC tracer package
subroutine ocmip2_cfc_end(CS)
  type(ocmip2_cfc_cs), pointer :: CS   !< The control structure returned by a
                                       !! previous call to register_OCMIP2_CFC.
!   This subroutine deallocates the memory owned by this module.
! Argument: CS - The control structure returned by a previous call to
!                register_OCMIP2_CFC.
  integer :: m

  if (associated(cs)) then
    if (associated(cs%CFC11)) deallocate(cs%CFC11)
    if (associated(cs%CFC12)) deallocate(cs%CFC12)

    deallocate(cs)
  endif
end subroutine ocmip2_cfc_end


!> \namespace mom_ocmip2_cfc
!!
!!   By Robert Hallberg, 2007
!!
!!     This module contains the code that is needed to set
!!   up and use CFC-11 and CFC-12 in a fully coupled or ice-ocean model
!!   context using the OCMIP2 protocols

end module mom_ocmip2_cfc