MOM_OCMIP2_CFC.F90

Go to the documentation of this file.

module mom_ocmip2_cfc
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, 2007                                           *
!*                                                                     *
!*    This file contains an example of the code that is needed to set  *
!*  up and use CFC-11 and CFC-12 in a fully coupled or ice-ocean model *
!*  context. There are 5 subroutines in this file.                     *
!*                                                                     *
!*    register_OCMIP2_CFC determines if the module is going to work,   *
!*  then makes several calls registering tracers to be advected and    *
!*  read from a restart file. it also sets various run-time parameters *
!*  for this module and sets up a "control structure" (CS) to store    *
!*  all information for this module.                                   *
!*                                                                     *
!*    initialize_OCMIP2_CFC initializes this modules arrays if they    *
!*  have not been found in a restart file.  It also determines which   *
!*  diagnostics will need to be calculated.                            *
!*                                                                     *
!*    OCMIP2_CFC_column_physics updates the CFC concentrations,        *
!*  applying everthing but horizontal advection and diffusion.         *
!*  Surface fluxes are applied inside an implicit vertical advection   *
!*  and diffusion tridiagonal solver, and any interior sources and     *
!*  sinks (not applicable for CFCs) would also be applied here.  This  *
!*  subroutine also sends out any requested interior diagnostics.      *
!*                                                                     *
!*    OCMIP2_CFC_surface_state calculates the information required     *
!*  from the ocean for the FMS coupler to calculate CFC fluxes.        *
!*                                                                     *
!*    OCMIP2_CFC_end deallocates the persistent run-time memory used   *
!*  by this module.                                                    *
!*                                                                     *
!*     A small fragment of the horizontal grid is shown below:         *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v, tr_ady, tr_dfy                        *
!*    j    x ^ x ^ x   At >:  u, tr_adx, tr_dfx                        *
!*    j    > o > o >   At o:  h, tr, CFC11, CFC12                      *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_ctrl
use mom_diag_to_z, only : register_z_tracer, diag_to_z_cs
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, read_data, slasher, vardesc, var_desc, query_vardesc
use mom_open_boundary, only : ocean_obc_type
use mom_restart, only : register_restart_field, query_initialized, mom_restart_cs
use mom_sponge, only : set_up_sponge_field, sponge_cs
use mom_time_manager, only : time_type, get_time
use mom_tracer_registry, only : register_tracer, tracer_registry_type
use mom_tracer_registry, only : add_tracer_diagnostics, add_tracer_obc_values
use mom_tracer_diabatic, only : tracer_vertdiff, applytracerboundaryfluxesinout
use mom_tracer_z_init, only : tracer_z_init
use mom_variables, only : surface
use mom_verticalgrid, only : verticalgrid_type

use coupler_util, only : extract_coupler_values, set_coupler_values
use coupler_util, only : ind_flux, ind_alpha, ind_csurf
use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux

implicit none ; private

#include <MOM_memory.h>

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


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

type p3d
  real, dimension(:,:,:), pointer :: p => null()
end type p3d

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 ! A pointer to the ocean model's clock.
  type(tracer_registry_type), pointer :: tr_reg => null()
  real, pointer, dimension(:,:,:) :: &
    cfc11 => null(), &     ! The CFC11 concentration in mol m-3.
    cfc12 => null(), &     ! The CFC12 concentration in mol m-3.
    cfc11_aux => null(), & ! The CFC11 and CFC12 concentrations, in mol m-3,
    cfc12_aux => null()    ! with values of thin layers masked out.
  ! In the following variables a suffix of _11 refers to CFC11 and _12 to CFC12.
  real :: a1_11, a2_11, a3_11, a4_11   ! Coefficients in the calculation of the
  real :: a1_12, a2_12, a3_12, a4_12   ! CFC11 and CFC12 Schmidt numbers, in
                                       ! units of ND, degC-1, degC-2, degC-3.
  real :: d1_11, d2_11, d3_11, d4_11   ! Coefficients in the calculation of the
  real :: d1_12, d2_12, d3_12, d4_12   ! CFC11 and CFC12 solubilities, in units
                                       ! of ND, K-1, log(K)^-1, K-2.
  real :: e1_11, e2_11, e3_11          ! More coefficients in the calculation of
  real :: e1_12, e2_12, e3_12          ! the CFC11 and CFC12 solubilities, in
                                       ! units of PSU-1, PSU-1 K-1, PSU-1 K-2.
  type(p3d), dimension(NTR) :: &
    tr_adx, &       ! Tracer zonal advective fluxes in mol s-1.
    tr_ady, &       ! Tracer meridional advective fluxes in mol s-1.
    tr_dfx, &       ! Tracer zonal diffusive fluxes in mol s-1.
    tr_dfy          ! Tracer meridional diffusive fluxes in mol s-1.
  real :: cfc11_ic_val = 0.0    ! The initial value assigned to CFC11.
  real :: cfc12_ic_val = 0.0    ! The initial value assigned to CFC12.
  real :: cfc11_land_val = -1.0 ! The values of CFC11 and CFC12 used where
  real :: cfc12_land_val = -1.0 ! land is masked out.
  logical :: mask_tracers  ! If true, tracers are masked out in massless layers.
  logical :: tracers_may_reinit  ! If true, tracers may go through the
                           ! initialization code if they are not found in the
                           ! restart files.
  character(len=16) :: cfc11_name, cfc12_name ! Variable names.

  integer :: ind_cfc_11_flux  ! Indices returned by aof_set_coupler_flux that
  integer :: ind_cfc_12_flux  ! are used to pack and unpack surface boundary
                              ! condition arrays.

  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                             ! timing of diagnostic output.
  type(mom_restart_cs), pointer :: restart_csp => null()
  integer :: id_cfc11, id_cfc12
  integer, dimension(NTR) :: id_tr_adx = -1, id_tr_ady = -1
  integer, dimension(NTR) :: id_tr_dfx = -1, id_tr_dfy = -1

  ! The following vardesc types contain a package of metadata about each tracer.
  type(vardesc) :: cfc11_desc, cfc12_desc
end type ocmip2_cfc_cs

contains

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.
! Arguments: HI - A horizontal index type structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
!  (in/out)  tr_Reg - A pointer to the tracer registry.
!  (in)      restart_CS - A pointer to the restart control structure.

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_OCMIP2_CFC" ! This module's name.
  character(len=200) :: inputdir ! The directory where NetCDF input files are.
  ! 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'
  real, dimension(:,:,:), pointer :: tr_ptr
  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.
  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)

  ! These calls obtain the indices for the CFC11 and CFC12 flux coupling.
  cs%ind_cfc_11_flux = 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")
  cs%ind_cfc_12_flux = 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")
  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 \n"//&
                 "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, "MASK_MASSLESS_TRACERS", cs%mask_tracers, &
                 "If true, the tracers are masked out in massless layer. \n"//&
                 "This can be a problem with time-averages.", default=.false.)
  call get_param(param_file, mdl, "TRACERS_MAY_REINIT", cs%tracers_may_reinit, &
                 "If true, tracers may go through the initialization code \n"//&
                 "if they are not found in the restart files.  Otherwise \n"//&
                 "it is a fatal error if tracers are not found in the \n"//&
                 "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)

  allocate(cs%CFC11(isd:ied,jsd:jed,nz)) ; cs%CFC11(:,:,:) = 0.0
  allocate(cs%CFC12(isd:ied,jsd:jed,nz)) ; cs%CFC12(:,:,:) = 0.0
  if (cs%mask_tracers) then
    allocate(cs%CFC11_aux(isd:ied,jsd:jed,nz)) ; cs%CFC11_aux(:,:,:) = 0.0
    allocate(cs%CFC12_aux(isd:ied,jsd:jed,nz)) ; cs%CFC12_aux(:,:,:) = 0.0
  endif

  ! 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 the restart file.
  call register_restart_field(tr_ptr, cs%CFC11_desc, &
                              .not.cs%tracers_may_reinit, restart_cs)
  ! Register CFC11 for horizontal advection & diffusion.
  call register_tracer(tr_ptr, cs%CFC11_desc, param_file, hi, gv, tr_reg, &
                       tr_desc_ptr=cs%CFC11_desc)
  ! Do the same for CFC12
  tr_ptr => cs%CFC12
  call register_restart_field(tr_ptr, cs%CFC12_desc, &
                              .not.cs%tracers_may_reinit, restart_cs)
  call register_tracer(tr_ptr, cs%CFC12_desc, param_file, hi, gv, tr_reg, &
                       tr_desc_ptr=cs%CFC12_desc)

  ! 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 NTR tracer fields in tr(:,:,:,:)
!! and it sets up the tracer output.
subroutine initialize_ocmip2_cfc(restart, day, G, GV, h, diag, OBC, CS, &
                                 sponge_CSp, diag_to_Z_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.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                  intent(in) :: h          !< Layer thicknesses, in H
                                                           !! (usually 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.
  type(diag_to_z_cs),             pointer    :: diag_to_Z_CSp !< A pointer to the control structure
                                                           !! for diagnostics in depth space.
!   This subroutine initializes the NTR tracer fields in tr(:,:,:,:)
! and it sets up the tracer output.

! Arguments: restart - .true. if the fields have already been read from
!                     a restart file.
!  (in)      day - Time of the start of the run.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      h - Layer thickness, in m or kg m-2.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (in)      OBC - This open boundary condition type specifies whether, where,
!                  and what open boundary conditions are used.
!  (in/out)  CS - The control structure returned by a previous call to
!                 register_OCMIP2_CFC.
!  (in/out)  sponge_CSp - A pointer to the control structure for the sponges, if
!                         they are in use.  Otherwise this may be unassociated.
!  (in/out)  diag_to_Z_Csp - A pointer to the control structure for diagnostics
!                            in depth space.
  logical :: from_file = .false.
  character(len=16) :: name     ! A variable's name in a NetCDF file.
  character(len=72) :: longname ! The long name of that variable.
  character(len=48) :: units    ! The dimensions of the variable.
  character(len=48) :: flux_units ! The units for tracer fluxes.
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
  integer :: IsdB, IedB, JsdB, JedB

  if (.not.associated(cs)) return
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  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, 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, cs)

  if (associated(obc)) then
  ! By default, all tracers have 0 concentration in their inflows. This may
  ! make the following calls are unnecessary.
  !  call add_tracer_OBC_values(trim(CS%CFC11_desc%name), CS%tr_Reg, 0.0)
  !  call add_tracer_OBC_values(trim(CS%CFC12_desc%name), CS%tr_Reg, 0.0)
  endif


  ! This needs to be changed if the units of tracer are changed above.
  if (gv%Boussinesq) then ; flux_units = "mol s-1"
  else ; flux_units = "mol m-3 kg s-1" ; endif

  do m=1,ntr
    ! Register the tracer advective and diffusive fluxes for potential
    ! diagnostic output.
    if (m==1) then
      ! Register CFC11 for potential diagnostic output.
      call query_vardesc(cs%CFC11_desc, name, units=units, longname=longname, &
                         caller="initialize_OCMIP2_CFC")
      cs%id_CFC11 = register_diag_field("ocean_model", trim(name), cs%diag%axesTL, &
          day, trim(longname) , trim(units))
      call register_z_tracer(cs%CFC11, trim(name), longname, units, &
                             day, g, diag_to_z_csp)
    elseif (m==2) then
      ! Register CFC12 for potential diagnostic output.
      call query_vardesc(cs%CFC12_desc, name, units=units, longname=longname, &
                         caller="initialize_OCMIP2_CFC")
      cs%id_CFC12 = register_diag_field("ocean_model", trim(name), cs%diag%axesTL, &
          day, trim(longname) , trim(units))
      call register_z_tracer(cs%CFC12, trim(name), longname, units, &
                             day, g, diag_to_z_csp)
    else
      call mom_error(fatal,"initialize_OCMIP2_CFC is only set up to work"//&
                           "with NTR <= 2.")
    endif

    cs%id_tr_adx(m) = register_diag_field("ocean_model", trim(name)//"_adx", &
        cs%diag%axesCuL, day, trim(longname)//" advective zonal flux" , &
        trim(flux_units))
    cs%id_tr_ady(m) = register_diag_field("ocean_model", trim(name)//"_ady", &
        cs%diag%axesCvL, day, trim(longname)//" advective meridional flux" , &
        trim(flux_units))
    cs%id_tr_dfx(m) = register_diag_field("ocean_model", trim(name)//"_dfx", &
        cs%diag%axesCuL, day, trim(longname)//" diffusive zonal flux" , &
        trim(flux_units))
    cs%id_tr_dfy(m) = register_diag_field("ocean_model", trim(name)//"_dfy", &
        cs%diag%axesCvL, day, trim(longname)//" diffusive zonal flux" , &
        trim(flux_units))
    if (cs%id_tr_adx(m) > 0) call safe_alloc_ptr(cs%tr_adx(m)%p,isdb,iedb,jsd,jed,nz)
    if (cs%id_tr_ady(m) > 0) call safe_alloc_ptr(cs%tr_ady(m)%p,isd,ied,jsdb,jedb,nz)
    if (cs%id_tr_dfx(m) > 0) call safe_alloc_ptr(cs%tr_dfx(m)%p,isdb,iedb,jsd,jed,nz)
    if (cs%id_tr_dfy(m) > 0) call safe_alloc_ptr(cs%tr_dfy(m)%p,isd,ied,jsdb,jedb,nz)

!    Register the tracer for horizontal advection & diffusion.
    if ((cs%id_tr_adx(m) > 0) .or. (cs%id_tr_ady(m) > 0) .or. &
        (cs%id_tr_dfx(m) > 0) .or. (cs%id_tr_dfy(m) > 0)) &
      call add_tracer_diagnostics(name, cs%tr_Reg, cs%tr_adx(m)%p, &
                                  cs%tr_ady(m)%p,cs%tr_dfx(m)%p,cs%tr_dfy(m)%p)
  enddo

end subroutine initialize_ocmip2_cfc
!>This subroutine initializes a tracer array.
subroutine init_tracer_cfc(h, tr, name, land_val, IC_val, G, CS)
  type(ocean_grid_type),                    intent(in)  :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)  :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: tr
  character(len=*),                         intent(in)  :: name
  real,                                     intent(in)  :: land_val, IC_val
  type(ocmip2_cfc_cs),                      pointer     :: CS

  ! 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)
      if (.not.ok) then
        ok = tracer_z_init(tr, h, cs%IC_file, trim(name), g)
        if (.not.ok) call mom_error(fatal,"initialize_OCMIP2_CFC: "//&
                "Unable to read "//trim(name)//" from "//&
                trim(cs%IC_file)//".")
      endif
    else
      call read_data(cs%IC_file, trim(name), tr, domain=g%Domain%mpp_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 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.
subroutine ocmip2_cfc_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, 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,
                                                           !! in m or kg m-2.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                      intent(in) :: h_new  !< Layer thickness after entrainment,
                                                           !! in 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, in 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, in m or kg m-2.
  type(forcing),                      intent(in) :: fluxes !< A structure containing pointers to any
                                                           !! possible forcing fields. Unused fields
                                                           !! have NULL ptrs.
  real,                               intent(in) :: dt     !< The amount of time covered by this
                                                           !! call, in s
  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
  real,                     optional,intent(in)  :: minimum_forcing_depth
!   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.

! Arguments: h_old -  Layer thickness before entrainment, in m or kg m-2.
!  (in)      h_new -  Layer thickness after entrainment, in m or kg m-2.
!  (in)      ea - an array to which the amount of fluid entrained
!                 from the layer above during this call will be
!                 added, in m or kg m-2.
!  (in)      eb - an array to which the amount of fluid entrained
!                 from the layer below during this call will be
!                 added, in m or kg m-2.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      dt - The amount of time covered by this call, in s.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 register_OCMIP2_CFC.
!
! 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]

  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, CFC12
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
  integer :: i, j, k, is, ie, js, je, nz, m

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  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 s-1].
  call extract_coupler_values(fluxes%tr_fluxes, cs%ind_cfc_11_flux, ind_flux, &
                              cfc11_flux, is, ie, js, je, -gv%Rho0)
  call extract_coupler_values(fluxes%tr_fluxes, cs%ind_cfc_12_flux, ind_flux, &
                              cfc12_flux, is, ie, js, je, -gv%Rho0)

  ! 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.
  if (cs%mask_tracers) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      if (h_new(i,j,k) < 1.1*gv%Angstrom) then
        cs%CFC11_aux(i,j,k) = cs%CFC11_land_val
        cs%CFC12_aux(i,j,k) = cs%CFC12_land_val
      else
        cs%CFC11_aux(i,j,k) = cfc11(i,j,k)
        cs%CFC12_aux(i,j,k) = cfc12(i,j,k)
      endif
    enddo ; enddo ; enddo
    if (cs%id_CFC11>0) call post_data(cs%id_CFC11, cs%CFC11_aux, cs%diag)
    if (cs%id_CFC12>0) call post_data(cs%id_CFC12, cs%CFC12_aux, cs%diag)
  else
    if (cs%id_CFC11>0) call post_data(cs%id_CFC11, cfc11, cs%diag)
    if (cs%id_CFC12>0) call post_data(cs%id_CFC12, cfc12, cs%diag)
  endif
  do m=1,ntr
    if (cs%id_tr_adx(m)>0) &
      call post_data(cs%id_tr_adx(m),cs%tr_adx(m)%p(:,:,:),cs%diag)
    if (cs%id_tr_ady(m)>0) &
      call post_data(cs%id_tr_ady(m),cs%tr_ady(m)%p(:,:,:),cs%diag)
    if (cs%id_tr_dfx(m)>0) &
      call post_data(cs%id_tr_dfx(m),cs%tr_dfx(m)%p(:,:,:),cs%diag)
    if (cs%id_tr_dfy(m)>0) &
      call post_data(cs%id_tr_dfy(m),cs%tr_dfy(m)%p(:,:,:),cs%diag)
  enddo

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, in H
                                                           !! (usually m or kg m-2).
  real, dimension(:),              intent(out)   :: stocks !< the mass-weighted integrated amount
                                                           !! of each tracer, in kg times
                                                           !! concentration units.
  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
! 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.

! Arguments: h - Layer thickness, in m or kg m-2.
!  (out)     stocks - the mass-weighted integrated amount of each tracer,
!                     in kg times concentration units.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 register_OCMIP2_CFC.
!  (out)     names - the names of the stocks calculated.
!  (out)     units - the units of the stocks calculated.
!  (in,opt)  stock_index - the coded index of a specific stock being sought.
! Return value: the number of stocks calculated here.

  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%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

subroutine ocmip2_cfc_surface_state(state, h, G, CS)
  type(ocean_grid_type),                    intent(in)    :: G    !< The ocean's grid structure.
  type(surface),                            intent(inout) :: state
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: h    !< Layer thicknesses, in H
                                                                  !! (usually m or kg m-2).
  type(ocmip2_cfc_cs),                      pointer       :: CS   !< The control structure returned
                                                                  !! by a previous call to
                                                                  !! register_OCMIP2_CFC.
!   This subroutine sets up the fields that the coupler needs to calculate the
! CFC fluxes between the ocean and atmosphere.
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (in)      h - Layer thickness, in m or kg m-2.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 register_OCMIP2_CFC.

  real, dimension(SZI_(G),SZJ_(G)) :: &
    CFC11_Csurf, &  ! The CFC-11 and CFC-12 surface concentrations times the
    CFC12_Csurf, &  ! Schmidt number term, both in mol m-3.
    CFC11_alpha, &  ! The CFC-11 solubility in mol m-3 pptv-1.
    CFC12_alpha     ! The CFC-12 solubility in mol m-3 pptv-1.
  real :: ta        ! Absolute sea surface temperature in units of dekaKelvin!?!
  real :: sal       ! Surface salinity in PSU.
  real :: SST       ! Sea surface temperature in degrees Celsius.
  real :: alpha_11  ! The solubility of CFC 11 in mol m-3 pptv-1.
  real :: alpha_12  ! The solubility of CFC 12 in 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, k, is, ie, js, je, m

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  if (.not.associated(cs)) return

  do j=js,je ; do i=is,ie
    ta = max(0.01, (state%SST(i,j) + 273.15) * 0.01) ! Why is this in hectoKelvin?
    sal = state%SSS(i,j) ; sst = 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 set_coupler_values(cfc11_alpha, state%tr_fields, cs%ind_cfc_11_flux, &
                          ind_alpha, is, ie, js, je)
  call set_coupler_values(cfc11_csurf, state%tr_fields, cs%ind_cfc_11_flux, &
                          ind_csurf, is, ie, js, je)
  call set_coupler_values(cfc12_alpha, state%tr_fields, cs%ind_cfc_12_flux, &
                          ind_alpha, is, ie, js, je)
  call set_coupler_values(cfc12_csurf, state%tr_fields, cs%ind_cfc_12_flux, &
                          ind_csurf, is, ie, js, je)

end subroutine ocmip2_cfc_surface_state

subroutine ocmip2_cfc_end(CS)
  type(ocmip2_cfc_cs), pointer :: CS
!   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)
    if (associated(cs%CFC11_aux)) deallocate(cs%CFC11_aux)
    if (associated(cs%CFC12_aux)) deallocate(cs%CFC12_aux)
    do m=1,ntr
      if (associated(cs%tr_adx(m)%p)) deallocate(cs%tr_adx(m)%p)
      if (associated(cs%tr_ady(m)%p)) deallocate(cs%tr_ady(m)%p)
      if (associated(cs%tr_dfx(m)%p)) deallocate(cs%tr_dfx(m)%p)
      if (associated(cs%tr_dfy(m)%p)) deallocate(cs%tr_dfy(m)%p)
    enddo

    deallocate(cs)
  endif
end subroutine ocmip2_cfc_end

end module mom_ocmip2_cfc