module field_ops
  !***************************************************************************
  !--- Perform operations on fields found in storage (interpolation, linear
  !--- combinations, units changes, etc) to prepare the set of variables that
  !--- will be transferred to either the ocean or the agcm.
  !--- This is done after all dependent fields have been received from other
  !--- component models and current versions of these fields are in storage.
  !---
  !--- Larry Solheim ...Sep,2015
  !***************************************************************************

  use com_cpl
  use ccc_esmf
  use ncdf
  use store
  use bulk
  use atm_f_data
  use timer
  use restart
  use spreading

  private

  public :: special_ops

  !--- Used for access to variable info (nx, ny, size, etc)
  type(cpl_vinfo_t), save :: curr_vinfo

  !--- Latent heat of sublimation for fresh ice or snow (J/kg)
  real(kind=8), parameter :: Lsub = 2.834e6

  !--- Conversion factor for CO2 in mass mixing ratio to CO2 in ppm
  real(kind=8), parameter :: co2_mmr2ppm = 1e6_8 * 28.97_8/44.0_8

  !--- The size of a 1D work array that will hold a single atm field
  integer, parameter :: sz_atm_data = (nlon_a + olap_a) * nlat_a

  !--- The size of a 1D work array that will hold a single ocn field
  integer, parameter :: sz_ocn_data = (nlon_o + olap_o) * nlat_o

  !--- Variables sent from the coupler to the agcm
  integer, save :: n_cpl2atm = 0
  character(32), pointer, save :: cpl2atm_name(:) => null()

  !--- Variables sent from the agcm to the coupler
  integer, save :: n_atm2cpl = 0
  character(32), pointer, save :: atm2cpl_name(:) => null()

  !--- Variables sent from the coupler to the ocean
  integer, save :: n_cpl2ocn = 0
  character(32), pointer, save :: cpl2ocn_name(:) => null()

  !--- Variables sent from the ocean to the coupler
  integer, save :: n_ocn2cpl = 0
  character(32), pointer, save :: ocn2cpl_name(:) => null()

  !--- Have these lists been initialized or not
  logical, save :: ops_initialized = .false.

  !--- ESMF related
  integer(ESMF_KIND_I) :: lnum
  integer(ESMF_KIND_I) :: rc

  integer, save :: nt_hist = 0

  contains

    subroutine read_snow_atm_from_file(snow)
      real(kind=8), pointer :: snow(:)

      !--- Local
      integer(4) :: fid, vid, ncstat, nwrds
      character(128) :: fname
      real(kind=8), pointer :: snow_tmp(:,:) => null()

      nullify(snow, snow_tmp)

      if ( store_GetIndex("SNOW_atm", other_gid) > 0 ) then
        !--- Read SNOW_atm from storage
        nwrds = store_readRecord("SNOW_atm", other_gid, snow, idx=1)

      else
        !--- Read SNOW_atm from and external file

        !--- This file contains a single monthly mean for SNOW_atm
        fname = "/fs/dev/crb/had_data/data/uncs_snow_atm_monthly.nc.001"
        call ncdf_open(fid, fname, "r")

        !--- Read SNOW_atm from the file
        call ncdf_read_var(fid, "SNOW_atm", snow_tmp)
        snow = reshape(snow_tmp, (/129*64/))

        call ncdf_close(fid)

        !--- Add this to storage
        call store_AppendRecord("SNOW_atm", other_gid, snow, create=.true.)
      endif

    end subroutine read_snow_atm_from_file

    !***************************************************************************
    !--- Initialize cpl2atm_name, cpl2ocn_name, atm2cpl_name, ocn2cpl_name
    !***************************************************************************
    subroutine field_ops_init()
      integer :: idx1, idx2, idx3, verbose=1
      character(32) :: tmp_names(200)
      logical :: found

      if (n_cpl2ocn == 0 .or. n_ocn2cpl == 0 .or. &
          n_cpl2atm == 0 .or. n_atm2cpl == 0) then
        !--- This initialization requires the event array to be defined
        if ( .not. associated(event) ) then
          write(6,*)"field_ops_init: event is not defined."
          call xit("field_ops_init",-1)
        endif
      endif

      if (n_cpl2ocn == 0) then
        !--- Define the list of variables sent from the coupler to the ocean

        tmp_names(:) = " "
        do idx1=1,nevents
          if ( event(idx1)%is_msg ) then
            !--- Ignore message variables
            cycle
          endif

          !--- Count the number of variables to be sent coupler to ocean
          !--- and populate tmp_names with their names
          if ( trim(adjustl(event(idx1)%get_model)) .eq. "cpl" .and. &
               trim(adjustl(event(idx1)%put_model)) .eq. "ocn" ) then
            if ( n_cpl2ocn == 0 ) then
              n_cpl2ocn = n_cpl2ocn + 1
              tmp_names(n_cpl2ocn) = trim(adjustl(event(idx1)%put_var))
            else
              !--- Make sure there are no duplicate names
              found = .false.
              do idx2=1,n_cpl2ocn
                if ( trim(adjustl(event(idx1)%put_var)) .eq. &
                     trim(adjustl(tmp_names(idx2))) ) then
                  !--- This is a duplicate variable name
                  found = .true.
                  exit
                endif
              enddo
              if ( .not. found ) then
                n_cpl2ocn = n_cpl2ocn + 1
                tmp_names(n_cpl2ocn) = trim(adjustl(event(idx1)%put_var))
              endif
            endif
          endif
        enddo

        if ( n_cpl2ocn > 0 ) then
          !--- Allocate and assign cpl2ocn_name
          allocate( cpl2ocn_name(n_cpl2ocn) )
          cpl2ocn_name(1:n_cpl2ocn) = tmp_names(1:n_cpl2ocn)
        else
          !--- Set n_cpl2ocn to -1 to indicate that there are no cpl to ocn
          !--- variables to be added to this list
          n_cpl2ocn = -1
        endif

        if ( verbose > 0 ) then
          write(6,'(a,i3)')'field_ops_init: n_cpl2ocn = ',n_cpl2ocn
          if ( n_cpl2ocn > 0 ) write(6,'(10a10)')cpl2ocn_name
          call flush(6)
        endif

      endif

      if (n_ocn2cpl == 0) then
        !--- Define the list of variables sent from the ocean to the coupler

        tmp_names(:) = " "
        do idx1=1,nevents
          if ( event(idx1)%is_msg ) then
            !--- Ignore message variables
            cycle
          endif

          !--- Count the number of variables to be sent ocean to coupler
          !--- and populate tmp_names with their names
          if ( trim(adjustl(event(idx1)%get_model)) .eq. "ocn" .and. &
               trim(adjustl(event(idx1)%put_model)) .eq. "cpl" ) then
            if ( n_ocn2cpl == 0 ) then
              n_ocn2cpl = n_ocn2cpl + 1
              tmp_names(n_ocn2cpl) = trim(adjustl(event(idx1)%put_var))
            else
              !--- Make sure there are no duplicate names
              found = .false.
              do idx2=1,n_ocn2cpl
                if ( trim(adjustl(event(idx1)%put_var)) .eq. &
                     trim(adjustl(tmp_names(idx2))) ) then
                  !--- This is a duplicate variable name
                  found = .true.
                  exit
                endif
              enddo
              if ( .not. found ) then
                n_ocn2cpl = n_ocn2cpl + 1
                tmp_names(n_ocn2cpl) = trim(adjustl(event(idx1)%put_var))
              endif
            endif
          endif
        enddo

        if ( n_ocn2cpl > 0 ) then
          !--- Allocate and assign ocn2cpl_name
          allocate( ocn2cpl_name(n_ocn2cpl) )
          ocn2cpl_name(1:n_ocn2cpl) = tmp_names(1:n_ocn2cpl)
        else
          !--- Set n_ocn2cpl to -1 to indicate that there are no ocn to cpl
          !--- variables to be added to this list
          n_ocn2cpl = -1
        endif

        if ( verbose > 0 ) then
          write(6,'(a,i3)')'field_ops_init: n_ocn2cpl = ',n_ocn2cpl
          if ( n_ocn2cpl > 0 ) write(6,'(10a10)')ocn2cpl_name
          call flush(6)
        endif

      endif

      if (n_cpl2atm == 0) then
        !--- Define the list of variables sent from the coupler to the agcm

        tmp_names(:) = " "
        do idx1=1,nevents
          if ( event(idx1)%is_msg ) then
            !--- Ignore message variables
            cycle
          endif

          !--- Count the number of variables to be sent coupler to atmosphere
          !--- and populate tmp_names with their names
          if ( trim(adjustl(event(idx1)%get_model)) .eq. "cpl" .and. &
               trim(adjustl(event(idx1)%put_model)) .eq. "atm" ) then
            if ( n_cpl2atm == 0 ) then
              n_cpl2atm = n_cpl2atm + 1
              tmp_names(n_cpl2atm) = trim(adjustl(event(idx1)%put_var))
            else
              !--- Make sure there are no duplicate names
              found = .false.
              do idx2=1,n_cpl2atm
                if ( trim(adjustl(event(idx1)%put_var)) .eq. &
                     trim(adjustl(tmp_names(idx2))) ) then
                  !--- This is a duplicate variable name
                  found = .true.
                  exit
                endif
              enddo
              if ( .not. found ) then
                n_cpl2atm = n_cpl2atm + 1
                tmp_names(n_cpl2atm) = trim(adjustl(event(idx1)%put_var))
              endif
            endif
          endif
        enddo

        if ( n_cpl2atm > 0 ) then
          !--- Allocate and assign cpl2atm_name
          allocate( cpl2atm_name(n_cpl2atm) )
          cpl2atm_name(1:n_cpl2atm) = tmp_names(1:n_cpl2atm)
        else
          !--- Set n_cpl2atm to -1 to indicate that there are no cpl to atm
          !--- variables to be added to this list
          n_cpl2atm = -1
        endif

        if ( verbose > 0 ) then
          write(6,'(a,i3)')'field_ops_init: n_cpl2atm = ',n_cpl2atm
          if ( n_cpl2atm > 0 ) write(6,'(10a10)')cpl2atm_name
          call flush(6)
        endif

      endif

      if (n_atm2cpl == 0) then
        !--- Define the list of variables sent from the agcm to the coupler

        tmp_names(:) = " "
        do idx1=1,nevents
          if ( event(idx1)%is_msg ) then
            !--- Ignore message variables
            cycle
          endif

          !--- Count the number of variables to be sent agcm to coupler
          !--- and populate tmp_names with their names
          if ( trim(adjustl(event(idx1)%get_model)) .eq. "atm" .and. &
               trim(adjustl(event(idx1)%put_model)) .eq. "cpl" ) then
            if ( n_atm2cpl == 0 ) then
              n_atm2cpl = n_atm2cpl + 1
              tmp_names(n_atm2cpl) = trim(adjustl(event(idx1)%put_var))
            else
              !--- Make sure there are no duplicate names
              found = .false.
              do idx2=1,n_atm2cpl
                if ( trim(adjustl(event(idx1)%put_var)) .eq. &
                     trim(adjustl(tmp_names(idx2))) ) then
                  !--- This is a duplicate variable name
                  found = .true.
                  exit
                endif
              enddo
              if ( .not. found ) then
                n_atm2cpl = n_atm2cpl + 1
                tmp_names(n_atm2cpl) = trim(adjustl(event(idx1)%put_var))
              endif
            endif
          endif
        enddo

        if ( n_atm2cpl > 0 ) then
          !--- Allocate and assign atm2cpl_name
          allocate( atm2cpl_name(n_atm2cpl) )
          atm2cpl_name(1:n_atm2cpl) = tmp_names(1:n_atm2cpl)
        else
          !--- Set n_atm2cpl to -1 to indicate that there are no agcm to cpl
          !--- variables to be added to this list
          n_atm2cpl = -1
        endif

        if ( verbose > 0 ) then
          write(6,'(a,i3)')'field_ops_init: n_atm2cpl = ',n_atm2cpl
          if ( n_atm2cpl > 0 ) write(6,'(10a10)')atm2cpl_name
          call flush(6)
        endif

      endif

      !--- Make it known that these lists have been initialized
      ops_initialized = .true.

    end subroutine field_ops_init

    !***************************************************************************
    !--- Handle non-transfer events
    !***************************************************************************
    function special_ops(event_id) result(found)

      implicit none

      !--- The event array index
      integer, intent(in) :: event_id
      logical :: found

      !--- Local
      integer :: idx, gid
      character(32) :: get_name, put_name
      character(len=256) :: strnga, strngb
      integer(impi) :: sz, tag, ierr
      logical :: comm_with_atm, comm_with_ocn

      found = .true.

      !--- Ensure the lists (cpl2atm_name, cpl2ocn_name, etc) containing variable
      !--- names of fields sent to the atm and to the ocn, are defined
      if ( .not. ops_initialized ) call field_ops_init()

      !--- Set local parameters from the event list
      get_name = trim(adjustl(event(event_id)%get_var))
      put_name = trim(adjustl(event(event_id)%put_var))

      comm_with_atm = .false.
      if ( mod(cpl_elapsed_time_secs, cpl_atm_freq) == 0 ) comm_with_atm = .true.
      comm_with_ocn = .false.
      if ( mod(cpl_elapsed_time_secs, cpl_ocn_freq) == 0 ) comm_with_ocn = .true.

      select case ( trim(get_name) )
        case ('start_cpl2atm')
          !--- The transfer from the coupler to the AGCM is about to start
          if ( comm_with_atm ) then
            !--- Stop timing ATM stop to ATM start
            call wall_time("ATM stop to ATM start", verbose=0)

            !--- Start timing the transfer from coupler to AGCM
            call wall_time("transfer CPL to ATM", reset=.true.)

            !--- Send a time and date string from the coupler to the atm
            !--- cpl_time_string is found in the com_cpl module
            call bcast_inter(cpl_time_string, cpl_master, "atm")
            !--- Send elasped coupler time in seconds from the coupler to the atm
            !--- cpl_elapsed_time_secs is found in the com_cpl module
            call bcast_inter(cpl_elapsed_time_secs, cpl_master, "atm")

            if ( atm_forcing_from_file ) then
              !--- Read atm forcing fields from a file
              call read_atm_forcing_data()
            endif

!DBG_bfb #define DBG_bfb
!DBG_bfb #ifdef DBG_bfb
!DBG_bfb strnga = " "
!DBG_bfb strngb = " "
!DBG_bfb write(strngb,'(i8.8)')cpl_elapsed_time_secs
!DBG_bfb if ( len_trim(env_runid) > 0 ) then
!DBG_bfb   strnga = "cplrs-before-cpl-to-atm_"//trim(strngb)//"_"//trim(env_runid)//".nc"
!DBG_bfb else
!DBG_bfb   strnga = "cplrs-before-cpl-to-atm_"//trim(strngb)//".nc"
!DBG_bfb endif
!DBG_bfb call restart_write(strnga, cpl_time_string)
!DBG_bfb #endif

          endif

        case ('stop_cpl2atm')
          !--- The transfer from the coupler to the AGCM has just completed
          if ( comm_with_atm .and. n_cpl2atm > 0 ) then
            !--- Stop timing the transfer from coupler to AGCM
            call wall_time("transfer CPL to ATM", verbose=0)

            !--- Remove the first field from the storage list associated with
            !--- each name for fields sent coupler to agcm
            !--- These fields will be associated with the ocn group because they
            !--- were interpolated using data from the ocean
            do idx=1,n_cpl2atm
              !--- Retain a maximum of 20 records in each list
              !--- Reset the group ID for certain variables
              gid = filter_group_ID(cpl2atm_name(idx), ocn_gid)
              call store_ShiftRecord( cpl2atm_name(idx), gid, max_size=2 )
            enddo

            !--- Start timing one AGCM coupling step
            call wall_time("ATM coupling step", reset=.true.)
          endif

        case ('start_cpl2ocn')
          !--- The transfer from the coupler to the ocean is about to start
          if ( comm_with_ocn ) then
            !--- Stop timing OCN stop to OCN start
            call wall_time("OCN stop to OCN start", verbose=0)

            !--- Start timing the transfer from coupler to OCN
            call wall_time("transfer CPL to OCN", reset=.true.)

            if ( .not. atm_forcing_from_file ) then
              !--- Send a time and date string from the coupler to the ocean
              !--- cpl_time_string is found in the com_cpl module
              call bcast_inter(cpl_time_string, cpl_master, "ocn")
              !--- Send elasped coupler time in seconds from the coupler to the ocean
              !--- cpl_elapsed_time_secs is found in the com_cpl module
              call bcast_inter(cpl_elapsed_time_secs, cpl_master, "ocn")
            endif

!DBG_bfb #ifdef DBG_bfb
!DBG_bfb strnga = " "
!DBG_bfb strngb = " "
!DBG_bfb write(strngb,'(i8.8)')cpl_elapsed_time_secs
!DBG_bfb if ( len_trim(env_runid) > 0 ) then
!DBG_bfb   strnga = "cplrs-before-cpl-to-ocn_"//trim(strngb)//"_"//trim(env_runid)//".nc"
!DBG_bfb else
!DBG_bfb   strnga = "cplrs-before-cpl-to-ocn_"//trim(strngb)//".nc"
!DBG_bfb endif
!DBG_bfb call restart_write(strnga, cpl_time_string)
!DBG_bfb #endif

          endif

        case ('stop_cpl2ocn')
          !--- The transfer from the coupler to the ocean has just completed
          if ( comm_with_ocn .and. n_cpl2ocn > 0 ) then
            !--- Stop timing the transfer from coupler to OCN
            call wall_time("transfer CPL to OCN", verbose=0)

            !--- Remove the first field from the storage list associated with
            !--- each name for fields sent coupler to ocean
            !--- These fields will be associated with the atm group because they
            !--- were interpolated using data from the agcm
            do idx=1,n_cpl2ocn
              !--- Retain a maximum of 20 records in each list
              gid = filter_group_ID(cpl2ocn_name(idx), atm_gid)
              call store_ShiftRecord( cpl2ocn_name(idx), gid, max_size=2 )
            enddo

            !--- Start timing one NEMO coupling step
            call wall_time("OCN coupling step", reset=.true.)
          endif

        case ('start_atm2cpl')
          !--- The transfer from the AGCM to the coupler is about to start
          if ( comm_with_atm ) then
            !--- Start timing the transfer from ATM to coupler
            call wall_time("transfer ATM to CPL", reset=.true.)

            !--- Start timing the transfer from ATM to coupler + prep ocn
            call wall_time("transfer ATM to CPL and prepare OCN fields", reset=.true.)
          endif

        case ('stop_atm2cpl')
          !--- The transfer from the AGCM to the coupler has just completed
          if ( comm_with_atm ) then

!DBG_bfb #ifdef DBG_bfb
!DBG_bfb strnga = " "
!DBG_bfb strngb = " "
!DBG_bfb write(strngb,'(i8.8)')cpl_elapsed_time_secs
!DBG_bfb if ( len_trim(env_runid) > 0 ) then
!DBG_bfb   strnga = "cplrs-after-atm-to-cpl_"//trim(strngb)//"_"//trim(env_runid)//".nc"
!DBG_bfb else
!DBG_bfb   strnga = "cplrs-after-atm-to-cpl_"//trim(strngb)//".nc"
!DBG_bfb endif
!DBG_bfb call restart_write(strnga, cpl_time_string)
!DBG_bfb #endif

            !--- Stop timing one AGCM coupling step
            call wall_time("ATM coupling step", verbose=0)

            !--- Stop timing the transfer from ATM to coupler
            call wall_time("transfer ATM to CPL", verbose=0)

            !--- Start timing prepare NEMO fields
            call wall_time("prepare NEMO fields", reset=.true.)

            !--- Define NEMO fields using the AGCM variables just transfered
            call prepare_nemo_fields()
            if ( n_atm2cpl > 0 ) then
              !--- Remove the first field from the storage list associated with
              !--- each name for fields sent agcm to coupler
              !--- These fields will be associated with the atm group because they
              !--- contain data received directly from the agcm
              do idx=1,n_atm2cpl
                !--- Retain a maximum of 2 records in each list
                gid = filter_group_ID(atm2cpl_name(idx), atm_gid)
                call store_ShiftRecord( atm2cpl_name(idx), gid, max_size=2 )
              enddo
            endif

            !--- Stop timing prepare NEMO fields
            call wall_time("prepare NEMO fields", verbose=0)

            !--- Stop timing the transfer from ATM to coupler + prep ocn
            call wall_time("transfer ATM to CPL and prepare OCN fields", verbose=0)

            !--- Start timing ATM stop to ATM start
            call wall_time("ATM stop to ATM start", reset=.true.)
          endif

        case ('start_ocn2cpl')
          !--- The transfer from the ocean to the coupler is about to start
          if ( comm_with_ocn ) then
            !--- Start timing the transfer from OCN to coupler
            call wall_time("transfer OCN to CPL", reset=.true.)

            !--- Start timing the transfer from OCN to coupler + prep atm
            call wall_time("transfer OCN to CPL and prepare ATM fields", reset=.true.)
          endif

        case ('stop_ocn2cpl')
          !--- The transfer from the ocean to the coupler has just completed
          if ( comm_with_ocn ) then

!DBG_bfb #ifdef DBG_bfb
!DBG_bfb strnga = " "
!DBG_bfb strngb = " "
!DBG_bfb write(strngb,'(i8.8)')cpl_elapsed_time_secs
!DBG_bfb if ( len_trim(env_runid) > 0 ) then
!DBG_bfb   strnga = "cplrs-after-ocn-to-cpl_"//trim(strngb)//"_"//trim(env_runid)//".nc"
!DBG_bfb else
!DBG_bfb   strnga = "cplrs-after-ocn-to-cpl_"//trim(strngb)//".nc"
!DBG_bfb endif
!DBG_bfb call restart_write(strnga, cpl_time_string)
!DBG_bfb #endif

            !--- Stop timing one NEMO coupling step
            call wall_time("OCN coupling step", verbose=0)

            !--- Stop timing the transfer from OCN to coupler
            call wall_time("transfer OCN to CPL", verbose=0)

            !--- Start timing prepare ATM fields
            call wall_time("prepare ATM fields", reset=.true.)

            !--- Define AGCM fields using the NEMO variables just transfered
            call prepare_agcm_fields()
            if ( n_ocn2cpl > 0 ) then
              !--- Remove the first field from the storage list associated with
              !--- each name for fields sent ocean to coupler
              !--- These fields will be associated with the ocn group because they
              !--- contain data received directly from the ocean
              do idx=1,n_ocn2cpl
                !--- Retain a maximum of 2 records in each list
                gid = filter_group_ID(ocn2cpl_name(idx), ocn_gid)
                call store_ShiftRecord( ocn2cpl_name(idx), gid, max_size=2 )
              enddo
            endif

            !--- Stop timing prepare ATM fields
            call wall_time("prepare ATM fields", verbose=0)

            !--- Stop timing the transfer from OCN to coupler + prep atm
            call wall_time("transfer OCN to CPL and prepare ATM fields", verbose=0)

            !--- Start timing OCN stop to OCN start
            call wall_time("OCN stop to OCN start", reset=.true.)
          endif

        case default
          found = .false.

      end select
    
    end function special_ops

    !***************************************************************************
    !--- Prepare fields for transfer from the coupler to NEMO
    !***************************************************************************
    subroutine prepare_nemo_fields()

      implicit none

      !--- Local
      integer :: idx1, idx2, ix, iy, nx, ny, nwrds, gid, grid_id, verbose=1
      real(kind=8) :: avalr8, bvalr8, cvalr8, rofo_val
      real(kind=8), pointer :: wrk1(:), wrk2(:), wrk3(:), wrk4(:), wrk5(:)
      character(64) :: src_fld_name
      integer(kind=8) :: store_time

      real(kind=8), dimension(:,:), pointer :: rnf_in, rnf_out, rnf_tmask, rnf_area

      !--- ESMF error string
      character(64) :: curr_method = "prepare_nemo_fields"

      !--- Grid area related
      type(ESMF_Field) :: area_field
      real(ESMF_KIND_R8), pointer :: area_data(:,:)

      !--- Should the CICE related flux fields (OTopMlt, OBotMlt, OIceEvap, OISubMSn)
      !--- passed to NEMO be grid box mean or per unit area
      logical :: cice_flux_per_unit_area

      nullify(wrk1, wrk2, wrk3, wrk4, wrk5)
      nullify(rnf_in, rnf_out, rnf_tmask, rnf_area, area_data)

      !--- This may go away in the future but keep it for the time being
      cice_flux_per_unit_area = .false.

      !--- The current time according to the coupler
      store_time = cpl_elapsed_time_secs

      !--- NOTE: In the following wrk1,2,3,4 will be reallocated by the
      !---       subroutines store_readRecord and regrid_Field (dst_name only)
      !---       on each call so it is not necessary to allocate them here

      !--- Add a storage field for each variable sent from the coupler to NEMO

      do idx1=1,n_cpl2ocn

        !--- This is the name associated with the ESMF field created in the ccc_esmf
        !--- module that will contain the source data
        src_fld_name = " "

        !--- The array wrk2 must contain the desired field on the NEMO grid
        !--- after the following case statement ends
        select case (trim(adjustl(cpl2ocn_name(idx1))))
          case ("O_OTaux1")
            !--- Atmosphere-Ocean wind stress (jpr_otx1)
            !--- Interpolated from AGCM variable OUFS_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("UFSO_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("UFSO_atm", -1)

            !--- Sign change for wind stress between ocean and atm
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("UFSO_atm", wrk1, "O_OTaux1", wrk2, windstress_remap)

          case ("O_OTauy1")
            !--- Atmosphere-Ocean wind stress (jpr_oty1)
            !--- Interpolated from AGCM variable OVFS_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("VFSO_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("VFSO_atm", -1)

            !--- Sign change for wind stress between ocean and atm
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("VFSO_atm", wrk1, "O_OTauy1", wrk2, windstress_remap)

          case ("O_OTaux2")
            !--- Atmosphere-Ocean wind stress (jpr_otx2)
            !--- Interpolated from AGCM variable OUFS_atm

            write(6,*)"prepare_nemo_fields: Invalid wind stress field request for O_OTaux2"
            write(6,*)"    If this is really what you want then edit field_ops.F90"
            call xit("prepare_nemo_fields",-1)

            !--- Get field data from internal storage
            nwrds = store_readRecord("OUFS_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("OUFS_atm", -1)

            !--- Sign change for wind stress between ocean and atm
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("OUFS_atm", wrk1, "O_OTaux2", wrk2, windstress_remap)

          case ("O_OTauy2")
            !--- Atmosphere-Ocean wind stress (jpr_oty2)
            !--- Interpolated from AGCM variable OVFS_atm

            write(6,*)"prepare_nemo_fields: Invalid wind stress field request for O_OTauy2"
            write(6,*)"    If this is really what you want then edit field_ops.F90"
            call xit("prepare_nemo_fields",-1)

            !--- Get field data from internal storage
            nwrds = store_readRecord("OVFS_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("OVFS_atm", -1)

            !--- Sign change for wind stress between ocean and atm
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("OVFS_atm", wrk1, "O_OTauy2", wrk2, windstress_remap)

          case ("O_ITaux1")
            !--- Atmosphere-Ice wind stress (jpr_itx1)
            !--- Interpolated from AGCM variable OUFS_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("UFSI_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("UFSI_atm", -1)

            !--- Sign change for wind stress between ocean and atm
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("UFSI_atm", wrk1, "O_ITaux1", wrk2, windstress_remap)

            if ( .false. .and. store_GetIndex("OIceFrc", ocn_gid) > 0 ) then
              !--- Get most recent sea ice fraction on the NEMO grid
              nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
              call  ReadRecord_nemo_err("OIceFrc", -1)

              !--- Determine grid box mean over ice covered cells
              wrk2 = wrk1*wrk2
            endif

          case ("O_ITauy1")
            !--- Atmosphere-Ice wind stress (jpr_ity1)
            !--- Interpolated from AGCM variable OVFS_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("VFSI_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("VFSI_atm", -1)

            !--- Sign change for wind stress between ocean and atm
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("VFSI_atm", wrk1, "O_ITauy1", wrk2, windstress_remap)

            if ( .false. .and. store_GetIndex("OIceFrc", ocn_gid) > 0 ) then
              !--- Get most recent sea ice fraction on the NEMO grid
              nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
              call  ReadRecord_nemo_err("OIceFrc", -1)

              !--- Determine grid box mean over ice covered cells
              wrk2 = wrk1*wrk2
            endif

          case ("O_ITaux2")
            !--- Atmosphere-Ice wind stress (jpr_itx2)
            !--- Interpolated from AGCM variable OUFS_atm

            write(6,*)"prepare_nemo_fields: Invalid wind stress field request for O_ITaux2"
            write(6,*)"    If this is really what you want then edit field_ops.F90"
            call xit("prepare_nemo_fields",-1)

            !--- Get field data from internal storage
            nwrds = store_readRecord("OUFS_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("OUFS_atm", -1)

            !--- Sign change for wind stress between ocean and atm
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("OUFS_atm", wrk1, "O_ITaux2", wrk2, windstress_remap)

            if ( .false. .and. store_GetIndex("OIceFrc", ocn_gid) > 0 ) then
              !--- Get most recent sea ice fraction on the NEMO grid
              nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
              call  ReadRecord_nemo_err("OIceFrc", -1)

              !--- Determine grid box mean over ice covered cells
              wrk2 = wrk1*wrk2
            endif

          case ("O_ITauy2")
            !--- Atmosphere-Ice wind stress (jpr_ity2)
            !--- Interpolated from AGCM variable OVFS_atm

            write(6,*)"prepare_nemo_fields: Invalid wind stress field request for O_ITauy2"
            write(6,*)"    If this is really what you want then edit field_ops.F90"
            call xit("prepare_nemo_fields",-1)

            !--- Get field data from internal storage
            nwrds = store_readRecord("OVFS_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("OVFS_atm", -1)

            !--- Sign change for wind stress between ocean and atm
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("OVFS_atm", wrk1, "O_ITauy2", wrk2, windstress_remap)

            if ( .false. .and. store_GetIndex("OIceFrc", ocn_gid) > 0 ) then
              !--- Get most recent sea ice fraction on the NEMO grid
              nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
              call  ReadRecord_nemo_err("OIceFrc", -1)

              !--- Determine grid box mean over ice covered cells
              wrk2 = wrk1*wrk2
            endif

          case ("O_QsrOce")
            !--- Solar heat flux over ocean (jpr_qsroce=13)
            !--- Interpolated from AGCM variable FSGO_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("FSGO_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("FSGO_atm", -1)

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("FSGO_atm", wrk1, "O_QsrOce", wrk2)

          case ("O_QsrIce")
            !--- Solar heat flux over ice (jpr_qsrice=14)
            !--- Interpolated from AGCM variable FSGI_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("FSGI_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("FSGI_atm", -1)

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("FSGI_atm", wrk1, "O_QsrIce", wrk2)

          case ("O_QsrMix")
            !--- Mixed Solar heat flux over ocean (jpr_qsrmix=15)
            !--- Interpolated from AGCM variable OFSG_atm

            !--- Get field data from internal storage
            !--- Solar over open water
            nwrds = store_readRecord("FSGO_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("FSGO_atm", -1)
            !--- Solar over ocean ice
            nwrds = store_readRecord("FSGI_atm", atm_gid, wrk2, idx=-1)
            call  ReadRecord_nemo_err("FSGI_atm", -1)
            !--- Sea ice fraction
            !--- SICN_atm was interpolated from ocean data at the previous time step
            !--- and so is associated with ocn_gid since that was the data origin
            !--- SICN is held constant in the AGCM over its coupling interval so this
            !--- value should be consistent with FSGO and FSGI
            nwrds = store_readRecord("SICN_atm", ocn_gid, wrk3, idx=-1)
            call  ReadRecord_nemo_err("SICN_atm", -1)

            !--- O_QsrMix = (1-SICN_atm)*FSGO_atm + SICN_atm*FSGI_atm
            wrk1 = (1.0_8-wrk3)*wrk1 + wrk3*wrk2

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("FSGO_atm", wrk1, "O_QsrMix", wrk2)

          case ("O_QnsOce")
            !--- Non-Solar heat flux over ocean (jpr_qnsoce=16)
            !--- Interpolated from AGCM variables (BEGO_atm - FSGO_atm)

            !--- Get field data from internal storage
            nwrds = store_readRecord("FSGO_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("FSGO_atm", -1)
            nwrds = store_readRecord("BEGO_atm", atm_gid, wrk2, idx=-1)
            call  ReadRecord_nemo_err("BEGO_atm", -1)

            src_fld_name = "BEGO_atm-FSGO_atm"

            !--- O_QnsOce = BEGO_atm - FSGO_atm
            wrk1 = wrk2 - wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("BEGO_atm", wrk1, "O_QnsOce", wrk2)

          case ("O_QnsIce")
            !--- Non-Solar heat flux over ice (jpr_qnsice=17)
            !--- Interpolated from AGCM variables (BEGI_atm - FSGI_atm)

            !--- Get field data from internal storage
            nwrds = store_readRecord("FSGI_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("FSGI_atm", -1)
            nwrds = store_readRecord("BEGI_atm", atm_gid, wrk2, idx=-1)
            call  ReadRecord_nemo_err("BEGI_atm", -1)

            src_fld_name = "BEGI_atm-FSGI_atm"

            !--- O_QnsIce = BEGI_atm - FSGI_atm
            wrk1 = wrk2 - wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("BEGI_atm", wrk1, "O_QnsIce", wrk2)

          case ("O_QnsMix")
            !--- Mixed Non-Solar heat flux over ocean (jpr_qnsmix=18)
            !--- Interpolated from AGCM variables (BEGO_atm - OFSG_atm)

            !--- Get field data from internal storage
            nwrds = store_readRecord("FSGO_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("FSGO_atm", -1)
            nwrds = store_readRecord("BEGO_atm", atm_gid, wrk2, idx=-1)
            call  ReadRecord_nemo_err("BEGO_atm", -1)

            nwrds = store_readRecord("FSGI_atm", atm_gid, wrk3, idx=-1)
            call  ReadRecord_nemo_err("FSGI_atm", -1)
            nwrds = store_readRecord("BEGI_atm", atm_gid, wrk4, idx=-1)
            call  ReadRecord_nemo_err("BEGI_atm", -1)

            !--- Sea ice fraction
            !--- SICN_atm was interpolated from ocean data at the previous time step
            !--- and so is associated with ocn_gid since that was the data origin
            !--- SICN is held constant in the AGCM over its coupling interval so this
            !--- value should be consistent with FSGO and FSGI
            nwrds = store_readRecord("SICN_atm", ocn_gid, wrk5, idx=-1)
            call  ReadRecord_nemo_err("SICN_atm", -1)

            src_fld_name = "O_QnsMix"

            !--- O_QnsOce = BEGO_atm - FSGO_atm
            !--- O_QnsIce = BEGI_atm - FSGI_atm
            !--- O_QnsMix = (1-SICN_atm)*O_QnsOce + SICN_atm*O_QnsIce
            wrk1 = wrk2 - wrk1
            wrk3 = wrk4 - wrk3
            wrk1 = (1.0_8-wrk5)*wrk1 + wrk5*wrk3

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("FSGO_atm", wrk1, "O_QnsMix", wrk2)

          case ("OTotRain")
            !--- Total liquid precipitation (jpr_rain=19)
            !--- Interpolated from AGCM variable RAIN_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("RAIN_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("RAIN_atm", -1)

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("RAIN_atm", wrk1, "OTotRain", wrk2)

            !--- Rain rate must be zero or positive
            !--- Interpolation can result in small negative values
            if ( any(wrk2 .lt. 0.0) ) then
              write(6,'(a)') &
                "prepare_nemo_fields: Adjusting negative interpolated rain values."
              where ( wrk2 < 0.0 ) wrk2 = 0.0
            endif

          case ("OTotSnow")
            !--- Total solid precipitation (jpr_snow=20)
            !--- Interpolated from AGCM variable ???? SNO_atm or SNOW_atm ????

            !--- Get field data from internal storage
            nwrds = store_readRecord("SNOW_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("SNOW_atm", -1)

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("SNOW_atm", wrk1, "OTotSnow", wrk2)

            !--- Snow rate must be zero or positive
            !--- Interpolation can result in small negative values
            if ( any(wrk2 .lt. 0.0) ) then
              write(6,'(a)') &
                "prepare_nemo_fields: Adjusting negative interpolated snow values."
              where ( wrk2 < 0.0 ) wrk2 = 0.0
            endif

          case ("OTotEvap")
            !--- Total evaporation (jpr_tevp=21)
            !--- Interpolated from AGCM variables (RAIN+SNOW-BWGO-BWGI)

            !--- Get field data from internal storage
            nwrds = store_readRecord("RAIN_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("RAIN_atm", -1)
            nwrds = store_readRecord("SNOW_atm", atm_gid, wrk2, idx=-1)
            call  ReadRecord_nemo_err("SNOW_atm", -1)
            nwrds = store_readRecord("BWGO_atm", atm_gid, wrk3, idx=-1)
            call  ReadRecord_nemo_err("BWGO_atm", -1)
            nwrds = store_readRecord("BWGI_atm", atm_gid, wrk4, idx=-1)
            call  ReadRecord_nemo_err("BWGI_atm", -1)
            nwrds = store_readRecord("SICN_atm", ocn_gid, wrk5, idx=-1)
            call  ReadRecord_nemo_err("SICN_atm", -1)

            src_fld_name = "RAIN_atm+SNOW_atm-BWGO_atm-BWGI_atm"

            !--- Evaluate total evaporation on the agcm grid
            !--- O_TotEvap = RAIN_atm + SNOW_atm - ( BWGI_atm*SICN_atm + BWGO_atm*(1-SICN_atm))
            wrk3 = wrk3*(1.0_8-wrk5)   !--- BWGO_atm*(1-SICN_atm)
            wrk4 = wrk4*wrk5           !--- BWGI_atm*SICN_atm
            wrk1 = wrk1 + wrk2 - wrk3 - wrk4

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("RAIN_atm", wrk1, "OTotEvap", wrk2)

          case ("OIceEvap")
            !--- Evaporation over ice (jpr_ievp=22)
            !--- Interpolated from AGCM variable (RAIN + SNOW - BWGI_atm)

            !--- Get field data from internal storage
            nwrds = store_readRecord("BWGI_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("BWGI_atm", -1)
            nwrds = store_readRecord("RAIN_atm", atm_gid, wrk2, idx=-1)
            call  ReadRecord_nemo_err("RAIN_atm", -1)
            nwrds = store_readRecord("SNOW_atm", atm_gid, wrk3, idx=-1)
            call  ReadRecord_nemo_err("SNOW_atm", -1)

            !--- RAIN + SNOW - BWGI
            wrk1 = wrk2 + wrk3 - wrk1

if (.false.) then
  !--- Previous definition
  nwrds = store_readRecord("HFLI_atm", atm_gid, wrk4, idx=-1)
  call  ReadRecord_nemo_err("HFLI_atm", -1)
  wrk4 = wrk4/Lsub
  if ( associated(wrk3) ) deallocate(wrk3)
  allocate(wrk3(size(wrk1)))
  wrk3 = wrk1 - wrk4
  avalr8 = minval(wrk3)
  bvalr8 = maxval(wrk3)
  where (wrk3 < 0.0_8) wrk3 = -wrk3
  where (wrk4 < 0.0_8) wrk4 = -wrk4
  cvalr8 = sum(wrk3)/sum(wrk4)
  write(6,'(a,3(a,g18.8))')trim(cpl_time_string), &
      "  OIceEvap: min(diff)=",avalr8,"  max(diff)=",bvalr8,"  relative(diff)=",cvalr8
  call flush(6)
endif

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("BWGI_atm", wrk1, "OIceEvap", wrk2)

            if ( cice_flux_per_unit_area ) then
              if ( store_GetIndex("OIceFrc", ocn_gid) > 0 ) then
                !--- Get sea ice fraction from the previous time step on the NEMO grid
                nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
                call  ReadRecord_nemo_err("OIceFrc", -1)

                !--- Multiply by ice fraction on the NEMO grid
                wrk2 = wrk1 * wrk2
              endif
            endif

          case ("OSubMPre")
            !--- Sublimation - liquid precip - solid precip  (jpr_sbpr=23)
            !--- Interpolated from AGCM variables (-BWGI_atm)

            !--- Get field data from internal storage
            nwrds = store_readRecord("BWGI_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("BWGI_atm", -1)

            src_fld_name = "-BWGI_atm"

            !--- Calculate submimation minus total precip over sea ice on the agcm grid as
            !--- -BWGI_atm
            wrk1 = -wrk1

if (.false.) then
  !--- Previous definition
  !--- BWGI = RAIN + SNOW - HFLI/Lsub = -OSubMPre
  nwrds = store_readRecord("HFLI_atm", atm_gid, wrk2, idx=-1)
  call  ReadRecord_nemo_err("HFLI_atm", -1)
  nwrds = store_readRecord("RAIN_atm", atm_gid, wrk3, idx=-1)
  call  ReadRecord_nemo_err("RAIN_atm", -1)
  nwrds = store_readRecord("SNOW_atm", atm_gid, wrk4, idx=-1)
  call  ReadRecord_nemo_err("SNOW_atm", -1)
  wrk4 = wrk2/Lsub - wrk3 - wrk4
  wrk3 = wrk1 - wrk4
  avalr8 = minval(wrk3)
  bvalr8 = maxval(wrk3)
  where (wrk3 < 0.0_8) wrk3 = -wrk3
  where (wrk4 < 0.0_8) wrk4 = -wrk4
  cvalr8 = sum(wrk3)/sum(wrk4)
  write(6,'(a,3(a,g18.8))')trim(cpl_time_string), &
      "  OSubMPre: min(diff)=",avalr8,"  max(diff)=",bvalr8,"  relative(diff)=",cvalr8
  call flush(6)
endif

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("BWGI_atm", wrk1, "OSubMPre", wrk2)

          case ("OISubMSn")
            !--- Ice freshwater budget  (jpr_semp=24)
            !--- Interpolated from AGCM variables (RAIN_atm - BWGI_atm)

            !--- Get field data from internal storage
            nwrds = store_readRecord("BWGI_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("BWGI_atm", -1)
            nwrds = store_readRecord( "RAIN_atm", atm_gid, wrk2, idx=-1)
            call  ReadRecord_nemo_err("RAIN_atm", -1)

            src_fld_name = "RAIN_atm - BWGI_atm"

            !--- Calculate ice freshwater budget on the agcm grid as
            !--- RAIN_atm - BWGI_atm
            wrk1 = wrk2 - wrk1

if (.false.) then
  !--- Previous definition
  nwrds = store_readRecord("HFLI_atm", atm_gid, wrk3, idx=-1)
  call  ReadRecord_nemo_err("HFLI_atm", -1)
  nwrds = store_readRecord("SNOW_atm", atm_gid, wrk4, idx=-1)
  call  ReadRecord_nemo_err("SNOW_atm", -1)
  wrk4 = wrk3/Lsub - wrk4
  wrk3 = wrk1 - wrk4
  avalr8 = minval(wrk3)
  bvalr8 = maxval(wrk3)
  where (wrk3 < 0.0_8) wrk3 = -wrk3
  where (wrk4 < 0.0_8) wrk4 = -wrk4
  cvalr8 = sum(wrk3)/sum(wrk4)
  write(6,'(a,3(a,g18.8))')trim(cpl_time_string), &
      "  OISubMSn: min(diff)=",avalr8,"  max(diff)=",bvalr8,"  relative(diff)=",cvalr8
  call flush(6)
endif

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("BWGI_atm", wrk1, "OISubMSn", wrk2)

            if ( cice_flux_per_unit_area ) then
              if ( store_GetIndex("OIceFrc", ocn_gid) > 0 ) then
                !--- Get sea ice fraction from the previous time step on the NEMO grid
                nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
                call  ReadRecord_nemo_err("OIceFrc", -1)

                !--- Multiply by ice fraction on the NEMO grid
                wrk2 = wrk1 * wrk2
              endif
            endif

          case ("OOEvaMPr")
            !--- Ocean freshwater budget  (jpr_oemp=25)
            !--- Interpolated from AGCM variable -BWGO_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("BWGO_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("BWGO_atm", -1)

            !--- BWGO is P-E whereas OOEvaMPr is E-P
            wrk1 = -wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("BWGO_atm", wrk1, "OOEvaMPr", wrk2)

          case ("O_Wind10")
            !--- 10 meter wind (jpr_w10m=26)
            !--- Interpolated from AGCM variable SWMX_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("SWMX_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("SWMX_atm", -1)

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("SWMX_atm", wrk1, "O_Wind10", wrk2)

          case ("O_dQnsdT")
            !--- non-Solar sensitivity over ice (jpr_dqnsdt=27)
            !--- Interpolated from AGCM variable SLIM_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("SLIM_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("SLIM_atm", -1)

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("SLIM_atm", wrk1, "O_dQnsdT", wrk2)

            if ( store_GetIndex("SLIMlw_atm", anygid) > 0 ) then
              !--- Also interpolate component terms of O_dQnsdT that are passed from the atm
              nwrds = store_readRecord("SLIMlw_atm", atm_gid, wrk1, idx=-1)
              if ( nwrds > 0 ) then
                call local_regrid_Field("SLIMlw_atm", wrk1, "O_dQlwdT", wrk3)
                gid = filter_group_ID("O_dQlwdT", atm_gid)
                call store_AppendRecord("O_dQlwdT", gid, wrk3, create=.true.)
              endif
            endif

            if ( store_GetIndex("SLIMsh_atm", anygid) > 0 ) then
              nwrds = store_readRecord("SLIMsh_atm", atm_gid, wrk1, idx=-1)
              if ( nwrds > 0 ) then
                call local_regrid_Field("SLIMsh_atm", wrk1, "O_dQshdT", wrk3)
                gid = filter_group_ID("O_dQshdT", atm_gid)
                call store_AppendRecord("O_dQshdT", gid, wrk3, create=.true.)
              endif
            endif

            if ( store_GetIndex("SLIMlh_atm", anygid) > 0 ) then
              nwrds = store_readRecord("SLIMlh_atm", atm_gid, wrk1, idx=-1)
              if ( nwrds > 0 ) then
                call local_regrid_Field("SLIMlh_atm", wrk1, "O_dQlhdT", wrk3)
                gid = filter_group_ID("O_dQlhdT", atm_gid)
                call store_AppendRecord("O_dQlhdT", gid, wrk3, create=.true.)
              endif
            endif

            !--- Also interpolate Sea Ice Temperature as known by the agcm, if present
            if ( store_GetIndex("GTICE_atm", anygid) > 0 ) then
              nwrds = store_readRecord("GTICE_atm", atm_gid, wrk1, idx=-1)
              if ( nwrds > 0 ) then
                if ( atm_forcing_from_file ) then
                  nwrds = store_readRecord("SICN_atm", anygid, wrk3, idx=-1)
                else
                  nwrds = store_readRecord("SICN_atm", ocn_gid, wrk3, idx=-1)
                endif
                if ( nwrds > 0 ) then
                  if ( .not. cpl_var_name_exists("GTICE_ocn") ) then
                    call define_cpl_var_by_name("GTICE_ocn", "ocn")
                  endif
                  !--- wrk1 = GTICE_atm * SICN_atm
                  wrk1 = wrk1*wrk3
                  call local_regrid_Field("GTICE_atm", wrk1, "GTICE_ocn", wrk4)
                  !--- wrk3 = SICN_atm
                  call local_regrid_Field("GTICE_atm", wrk3, "GTICE_ocn", wrk5)
                  where (wrk5 >  0.0) wrk4 = wrk4/wrk5
                  where (wrk5 <= 0.0) wrk4 = 0.0_8
                  gid = filter_group_ID("GTICE_ocn", atm_gid)
                  call store_AppendRecord("GTICE_ocn", gid, wrk4, create=.true.)
                endif
              endif
            endif

            if ( bulk_in_cpl ) then
              !--- Evaluate O_dQnsdT using a bulk formulation
              !--- Overwrite wrk2 with the value from bulk_dqnsdt
              write(6,*)"prepare_nemo_fields: Using bulk formulation for O_dQnsdT"
              call flush(6)
              call bulk_dqnsdt(wrk2, wrk3, wrk4, wrk5)

              if ( .not. cpl_var_name_exists("bulk_dQnsdT") ) then
                call define_cpl_var_by_name("bulk_dQnsdT", "ocn")
              endif
              gid = filter_group_ID("bulk_dQnsdT", atm_gid)
              call store_AppendRecord("bulk_dQnsdT", gid, wrk2, create=.true.)

              if ( .not. cpl_var_name_exists("bulk_dQlhdT") ) then
                call define_cpl_var_by_name("bulk_dQlhdT", "ocn")
              endif
              gid = filter_group_ID("bulk_dQlhdT", atm_gid)
              call store_AppendRecord("bulk_dQlhdT", gid, wrk3, create=.true.)

              if ( .not. cpl_var_name_exists("bulk_dQlwdT") ) then
                call define_cpl_var_by_name("bulk_dQlwdT", "ocn")
              endif
              gid = filter_group_ID("bulk_dQlwdT", atm_gid)
              call store_AppendRecord("bulk_dQlwdT", gid, wrk4, create=.true.)

              if ( .not. cpl_var_name_exists("bulk_dQshdT") ) then
                call define_cpl_var_by_name("bulk_dQshdT", "ocn")
              endif
              gid = filter_group_ID("bulk_dQshdT", atm_gid)
              call store_AppendRecord("bulk_dQshdT", gid, wrk5, create=.true.)
            endif

          case ("O_Runoff")
            !--- Runoff (jpr_rnf=28)

            !--- Get field data from internal storage
            nwrds = store_readRecord("RIVO_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("RIVO_atm", -1)

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("RIVO_atm", wrk1, "O_Runoff", wrk2)

            if ( size(wrk2) /= nemo_jpiglo*nemo_jpjglo ) then
              write(6,*)"prepare_nemo_fields: runoff_spreading grid size mismatch"
              write(6,*)"    expected nx*ny = ",nemo_jpiglo*nemo_jpjglo
              write(6,*)"       found nx*ny = ",size(wrk2)
              call xit("prepare_nemo_fields",-2)
            endif

            if ( associated(rnf_in) ) deallocate(rnf_in)
            allocate( rnf_in(nemo_jpiglo,nemo_jpjglo) )
            if ( associated(rnf_out) ) deallocate(rnf_out)
            allocate( rnf_out(nemo_jpiglo,nemo_jpjglo) )
            if ( associated(rnf_tmask) ) deallocate(rnf_tmask)
            allocate( rnf_tmask(nemo_jpiglo,nemo_jpjglo) )
            if ( associated(rnf_area) ) deallocate(rnf_area)
            allocate( rnf_area(nemo_jpiglo,nemo_jpjglo) )

            rnf_tmask = nemo_tmask
            rnf_area  = nemo_e1t * nemo_e2t

write(6,*)"O_Runoff: Before runoff_spreading min=",minval(wrk2),"  max=",maxval(wrk2)

            rnf_in = reshape(wrk2, (/nemo_jpiglo,nemo_jpjglo/) )
            call runoff_spreading(rnf_in, rnf_tmask, rnf_area, rnf_out)
            wrk2 = reshape(rnf_out, (/nemo_jpiglo*nemo_jpjglo/) )

write(6,*)"O_Runoff: After runoff_spreading min=",minval(wrk2),"  max=",maxval(wrk2)
call flush(6)

            !--- Clean up
            deallocate( rnf_in, rnf_out, rnf_tmask, rnf_area )

            if ( zero_runoff_sent_to_ocean ) then
              !--- The user has requested that zero runoff sent to the ocean
              wrk2(:) = 0.0
              write(6,*)"prepare_nemo_fields: O_Runoff is zeroed before being sent to NEMO"
              write(6,*)"prepare_nemo_fields: Avg value of RIVO_atm sent from the AGCM is ",sum(wrk1)/size(wrk1)
              call flush(6)
            endif

          case ("OCalving")
            !--- Calving (jpr_cal=29)
            !--- Calving is counted in runoff so always send zero

            if ( associated(wrk2) ) deallocate(wrk2)
            allocate( wrk2(sz_ocn_data) )
            wrk2(:) = 0.0

          case ("O_TauMod")
            !--- Wind stress amplitude (jpr_taum=30)
            !--- Interpolated from AGCM variables sqrt(UFSO_atm^2 + VFSO_atm^2)

            !--- Get field data from internal storage
            if ( store_GetIndex("UFSO_atm", atm_gid) > 0 ) then
              !--- UFSO_atm exists in storage, assume VFSO_atm does as well
              nwrds = store_readRecord("UFSO_atm", atm_gid, wrk1, idx=-1)
              call  ReadRecord_nemo_err("UFSO_atm", -1)
              nwrds = store_readRecord("VFSO_atm", atm_gid, wrk2, idx=-1)
              call  ReadRecord_nemo_err("VFSO_atm", -1)

              !--- Calculate wind stress amplitude on the agcm grid
#undef taumod_with_ice
#ifdef taumod_with_ice
              if ( store_GetIndex("UFSI_atm", atm_gid) > 0 ) then
                !--- UFSI_atm exists in storage, assume VFSI_atm does as well
                nwrds = store_readRecord("UFSI_atm", atm_gid, wrk3, idx=-1, pfx=" ")
                call  ReadRecord_nemo_err("UFSI_atm", -1)
                nwrds = store_readRecord("VFSI_atm", atm_gid, wrk4, idx=-1, pfx=" ")
                call  ReadRecord_nemo_err("VFSI_atm", -1)
                !--- This is SICN interpolated from ocean data at the previous time step
                nwrds = store_readRecord("SICN_atm", ocn_gid, wrk5, idx=-1, pfx=" ")
                call  ReadRecord_nemo_err("SICN_atm", -1)
                wrk1 = (1.0-wrk5) * sqrt(wrk1**2 + wrk2**2) + &
                            wrk5  * sqrt(wrk3**2 + wrk4**2)
                src_fld_name = "Mixed sqrt(UFSOI_atm^2+VFSOI_atm^2)"
              else
                wrk1 = sqrt(wrk1**2 + wrk2**2)
                src_fld_name = "sqrt(UFSO_atm^2+VFSO_atm^2)"
              endif
#else
              wrk1 = sqrt(wrk1**2 + wrk2**2)
              src_fld_name = "sqrt(UFSO_atm^2+VFSO_atm^2)"
#endif

            else if ( store_GetIndex("OUFS_atm", atm_gid) > 0 ) then
              !--- OUFS_atm and OVFS_atm are non-tiled fields which were use before
              !--- UFSO_atm and VFSO_atm (tiled fields) were avaialable

              !--- OUFS_atm exists in storage, assume OVFS_atm does as well
              nwrds = store_readRecord("OUFS_atm", atm_gid, wrk1, idx=-1)
              call  ReadRecord_nemo_err("OUFS_atm", -1)
              nwrds = store_readRecord("OVFS_atm", atm_gid, wrk2, idx=-1)
              call  ReadRecord_nemo_err("OVFS_atm", -1)

              !--- Calculate wind stress amplitude on the agcm grid
              wrk1 = sqrt(wrk1**2 + wrk2**2)

              src_fld_name = "sqrt(OUFS_atm^2+OVFS_atm^2)"

            else
              write(6,*)"prepare_nemo_fields: Unable to process O_TauMod."
              write(6,*)"prepare_nemo_fields: Neither UFSO_atm nor OUFS_atm exist in storage."
              call xit("prepare_nemo_fields",-2)
            endif

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("UFSO_atm", wrk1, "O_TauMod", wrk2)

          case ("O_AtmCO2")
            !--- CO2 concentration sent from the atm to NEMO (jpr_co2=31)

            !--- Get field data from internal storage
            nwrds = store_readRecord("CO2_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("CO2_atm", -1)

write(6,*)"CO2_atm: Received from AGCM  min=",minval(wrk1),"  max=",maxval(wrk1)

            !--- Convert from mass mixing ratio to ppm
            wrk1 = co2_mmr2ppm * wrk1

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("CO2_atm", wrk1, "O_AtmCO2", wrk2)

write(6,*)"O_AtmCO2: Sent to NEMO  min=",minval(wrk2),"  max=",maxval(wrk2)
call flush(6)

          case ("O_MSLP")
            !--- Mean sea level pressure sent from the atm to NEMO (jpr_mslp=34)

            !--- Get field data from internal storage
            nwrds = store_readRecord("PMSL_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("PMSL_atm", -1)

write(6,*)"PMSL_atm: Received from AGCM  min=",minval(wrk1),"  max=",maxval(wrk1)

            !--- Convert from hPa to atmospheres
            wrk1 = wrk1/1013.25_8

            !--- Interpolate this onto the OCN grid
            call local_regrid_Field("PMSL_atm", wrk1, "O_MSLP", wrk2)

write(6,*)"O_MSLP: Sent to NEMO  min=",minval(wrk2),"  max=",maxval(wrk2)
call flush(6)

          case ("OTopMlt")
            !--- top melt (jpr_topm=32)
            !--- net flux minus conductive flux at the top of the ice
            !--- Interpolated from AGCM variable BEGI_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("BEGI_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("BEGI_atm", -1)

            !--- Interpolate this onto the NEMO grid
            call local_regrid_Field("BEGI_atm", wrk1, "OTopMlt", wrk2)

            if ( cice_flux_per_unit_area ) then
              if ( store_GetIndex("OIceFrc", ocn_gid) > 0 ) then
                !--- Get sea ice fraction from the previous time step on the NEMO grid
                nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
                call  ReadRecord_nemo_err("OIceFrc", -1)

                !--- Multiply by ice fraction on the NEMO grid
                wrk2 = wrk1 * wrk2
              endif
            endif

          case ("OBotMlt")
            !--- bottom melt (jpr_botm=33)
            !--- conductive flux at the top of the ice (positive down)
            !--- Interpolated from AGCM variable -HSEA_atm

            !--- Get field data from internal storage
            nwrds = store_readRecord("HSEA_atm", atm_gid, wrk1, idx=-1)
            call  ReadRecord_nemo_err("HSEA_atm", -1)

            !??? Should this be -HSEA_atm ???
            src_fld_name = "HSEA_atm"

            !--- Interpolate this onto the NEMO grid
            call local_regrid_Field("HSEA_atm", wrk1, "OBotMlt", wrk2)

            if ( cice_flux_per_unit_area ) then
              if ( store_GetIndex("OIceFrc", ocn_gid) > 0 ) then
                !--- Get sea ice fraction from the previous time step on the NEMO grid
                nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
                call  ReadRecord_nemo_err("OIceFrc", -1)

                !--- Multiply by ice fraction on the NEMO grid
                wrk2 = wrk1 * wrk2
              endif
            endif

          case default
            !--- Only the field names listed above are currently recognized
            !--- Add additional cases as required
            write(6,'(3a)') &
              'prepare_nemo_fields: Unknown field name --> ',trim(cpl2ocn_name(idx1)),' <--'
            call flush(6)
            call xit("prepare_nemo_fields",-1)
        end select

        if ( verbose > 5 ) then
          call print_var_stats(wrk2, name=cpl2ocn_name(idx1), &
                               pfx="prepare_nemo_fields: "//trim(cpl2ocn_name(idx1)))
        endif

        !--- Add this data to internal storage
        !--- These fields are associated with the atm group because they were
        !--- interpolated from data that was received from the agcm
        gid = filter_group_ID(cpl2ocn_name(idx1), atm_gid)
        call store_AppendRecord(cpl2ocn_name(idx1), gid, wrk2, create=.true.)

      enddo

      !--- Clean up
      if ( associated(wrk1) ) deallocate(wrk1)
      if ( associated(wrk2) ) deallocate(wrk2)
      if ( associated(wrk3) ) deallocate(wrk3)
      if ( associated(wrk4) ) deallocate(wrk4)
      if ( associated(wrk5) ) deallocate(wrk5)

      contains

        subroutine local_regrid_field(src_name, src_data, dst_name, dst_data, methin)
          character(*) :: src_name, dst_name
          character(*), optional :: methin
          character(len=32) :: rmmethod
          real(kind=8), intent(inout), pointer :: src_data(:), dst_data(:)

          if ( len_trim(src_fld_name) < 1 ) then
            src_fld_name = trim(src_name)
          endif

          if (present(methin)) then
             rmmethod=trim(methin)
          else
             rmmethod="conservative"
          endif


          !--- Conservative-fracarea remap from AGCM to NEMO grid with land values on
          !--- the NEMO grid masked (land does not take part in the remap calculation)
          call regrid_Field(src_name, src_data, dst_name, dst_data, &
                            method=trim(rmmethod), norm_type="fracarea", &
                            src_fld_name=trim(src_fld_name), &
                            dst_mask_values=(/0/), missing_value=0.0_8)

        end subroutine local_regrid_field

        subroutine ReadRecord_nemo_err(name, err)
          character(*) :: name
          integer :: err
          if ( nwrds < 0 ) then
            !--- The field was not found in storage
            write(6,*)"prepare_nemo_fields: ",trim(cpl2ocn_name(idx1)), &
                " ...",trim(name)," is missing from storage."
            call xit("prepare_nemo_fields",err)
          endif
        end subroutine ReadRecord_nemo_err

    end subroutine prepare_nemo_fields

    !***************************************************************************
    !--- Prepare fields for transfer from the coupler to the AGCM
    !***************************************************************************
    subroutine prepare_agcm_fields()

      implicit none

      !--- Local
      integer :: idx1, nwrds, gid, verbose=1
      integer :: nx_hist, ny_hist
      integer(kind=8) :: store_time
      real(kind=8), pointer :: wrk1(:), wrk2(:), wrk3(:), wrk4(:)
      real(kind=8) :: sea_ice_density, snow_density
      character(64) :: src_fld_name
      character(256) :: strnga, strngb

      nullify(wrk1, wrk2, wrk3, wrk4)

      store_time = cpl_elapsed_time_secs

      !--- NOTE: In the following wrk1,2 will be reallocated by the
      !---       subroutines store_readRecord and regrid_Field (dst_name only)
      !---       on each call so it is not necessary to allocate them here

      !--- Add a storage field for each variable sent from cpl to atm

      do idx1=1,n_cpl2atm

        !--- This is the name associated with the ESMF field created in the ccc_esmf
        !--- module that will contain the source data
        src_fld_name = " "

        select case (trim(adjustl(cpl2atm_name(idx1))))
          case ("GT_atm")
            !--- Ground temperature
            !--- Interpolated from NEMO variable O_SSTSST or O_TepMix

            !--- nemo_namsbc_cpl_cldes(1) is sn_snd_temp%cldes sent from NEMO
            select case ( trim(adjustl(nemo_namsbc_cpl_cldes(1))) )
              case ("oce only")
                !--- NEMO sends O_SSTSST = tsn+273.15
                !--- O_SSTSST is SST defined over the entire ocean including under sea ice
                !xxx nwrds = store_readRecord("O_SSTSST", ocn_gid, wrk1, time=store_time)
                nwrds = store_readRecord("O_SSTSST", ocn_gid, wrk1, idx=-1)
                call  ReadRecord_agcm_err("O_SSTSST", -1)
                src_fld_name = "O_SSTSST"

              case ("weighted oce and ice")
                !--- NEMO sends O_SSTSST = (tsn+273.15)*(1-fr_i) and O_TepIce = tn_ice*fr_i
                !--- Combine these 2 fields to get a mixed surface temperature which is
                !--- sea surface temperature over open water and ice surface temperature
                !--- over ice covered areas
                !xxx nwrds = store_readRecord("O_SSTSST", ocn_gid, wrk1, time=store_time)
                nwrds = store_readRecord("O_SSTSST", ocn_gid, wrk1, idx=-1)
                call  ReadRecord_agcm_err("O_SSTSST", -1)
                !xxx nwrds = store_readRecord("O_TepIce", ocn_gid, wrk2, time=store_time)
                nwrds = store_readRecord("O_TepIce", ocn_gid, wrk2, idx=-1)
                call  ReadRecord_agcm_err("O_TepIce", -1)
                wrk1 = wrk1 + wrk2
                src_fld_name = "O_SSTSST and O_TepIce"

              case ("oce and ice")
                !--- NEMO sends O_SSTSST = (tsn+273.15) and O_TepIce = tn_ice
                !--- Use O_SSTSST as is
                !xxx nwrds = store_readRecord("O_SSTSST", ocn_gid, wrk1, time=store_time)
                nwrds = store_readRecord("O_SSTSST", ocn_gid, wrk1, idx=-1)
                call  ReadRecord_agcm_err("O_SSTSST", -1)
                src_fld_name = "O_SSTSST"

              case ("mixed oce-ice")
                !--- NEMO sends O_TepMix = (tsn+273.15)*(1-fr_i) + tn_ice*fr_i
                !--- This is a mixed surface temperature which is sea surface
                !--- temperature over open water and ice surface temperature
                !--- over ice covered areas
                !xxx nwrds = store_readRecord("O_TepMix", ocn_gid, wrk1, time=store_time)
                nwrds = store_readRecord("O_TepMix", ocn_gid, wrk1, idx=-1)
                call  ReadRecord_agcm_err("O_TepMix", -1)
                src_fld_name = "O_TepMix"

              case default
                write(6,*)"prepare_agcm_fields: Invalid description for O_SSTSST --> ", &
                          trim(nemo_namsbc_cpl_cldes(1)),' <--'
                call flush(6)
                call xit("prepare_agcm_fields",-1)

            end select


#undef DBG_bfb
#ifdef DBG_bfb
call dump_farray("O_SSTSST_before_regrid",wrk1,idx=cpl_elapsed_time_secs)
#endif

            !--- Interpolate this onto the ATM grid
            call local_regrid_Field("O_SSTSST", wrk1, "GT_atm", wrk2)


#ifdef DBG_bfb
call dump_farray("GT_atm_after_regrid",wrk2,idx=cpl_elapsed_time_secs)
#endif

!???            !--- Post process GT_atm after interpolation in certain cases
!???            select case ( trim(adjustl(nemo_namsbc_cpl_cldes(1))) )
!???              case ("weighted oce and ice", "mixed oce-ice")
!???                !--- GT_atm will contain a mixed SST and sea ice temperature
!???                !--- Untangle this so that GT_atm contains only SST
!???!---TODO---
!???!--- Maybe it is better to do before interpolation in the case statement above
!???                write(6,*)"prepare_agcm_fields: Invalid description for O_SSTSST --> ", &
!???                          trim(nemo_namsbc_cpl_cldes(1)),' <--  Not yet implimented.'
!???                call flush(6)
!???                call xit("prepare_agcm_fields",-1)
!???
!???                !--- Get the value of GT_atm from the previous time step
!???                nwrds = store_readRecord("GT_atm", atm_gid, wrk1, &
!???                                         time=store_time-cpl_time_step_secs)
!???                call  ReadRecord_agcm_err("GT_atm", -1)
!???            end select

!xxxCanOM4          case ("GC_atm")
!xxxCanOM4            !--- Ground cover mask
!xxxCanOM4            !--- This is not used when coupled with NEMO
!xxxCanOM4
!xxxCanOM4            !--- Get the most recent GC_atm from internal storage
!xxxCanOM4            !--- No interpolation is done on GC_atm
!xxxCanOM4            nwrds = store_readRecord("GC_atm", atm_gid, wrk2, idx=-1)
!xxxCanOM4            call  ReadRecord_agcm_err("GC_atm", -1)

          case ("SICN_atm")
            !--- Sea ice fraction
            !--- Interpolated from NEMO variable OIceFrc

            if ( nemo_nn_ice > 0 ) then
              !--- Sea ice is used in NEMO and should be found in storage
              !--- Get field data from internal storage
              !xxx nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, time=store_time)
              nwrds = store_readRecord("OIceFrc", ocn_gid, wrk1, idx=-1)
              call  ReadRecord_agcm_err("OIceFrc", -1)

#ifdef DBG_bfb
call dump_farray("OIceFrc_before_regrid",wrk1,idx=cpl_elapsed_time_secs)
#endif

              !--- Interpolate this onto the ATM grid
              call local_regrid_Field("OIceFrc", wrk1, "SICN_atm", wrk2)

#ifdef DBG_bfb
call dump_farray("SICN_atm_after_regrid",wrk2,idx=cpl_elapsed_time_secs)
#endif

            else
              !--- Sea ice is not used in NEMO
              !--- Create an array that is identically zero for SICN_atm
              if ( associated(wrk2) ) deallocate(wrk2)
              allocate( wrk2(sz_atm_data) )
              wrk2(:) = 0.0
            endif

            !--- Interpolation can lead to small negative values
            where ( wrk2 < 0.0 ) wrk2 = 0.0

          case ("SIC_atm")
            !--- Interpolated from NEMO variable OIceTck
            !--- if sn_snd_thick%cldes = 'ice and snow' then
            !---     OIceTck = Sea ice thickness
            !--- if sn_snd_thick%cldes = 'weighted ice and snow' then
            !---     OIceTck = Cell average sea ice thickness
            !--- if sn_snd_thick%cldes = 'weighted iwe and swe'
            !---     OIceTck = Sea ice water equivalent

            !--- Get field data from internal storage
            !xxx nwrds = store_readRecord("OIceTck", ocn_gid, wrk1, time=store_time)
            nwrds = store_readRecord("OIceTck", ocn_gid, wrk1, idx=-1)
            call  ReadRecord_agcm_err("OIceTck", -1)

            !--- nemo_namsbc_cpl_cldes(3) is sn_snd_thick%cldes sent from NEMO
            select case ( trim(adjustl(nemo_namsbc_cpl_cldes(3))) )
              case ("ice and snow")
                !--- OIceTck is the non-averaged sea ice thickness (ht_i in NEMO)
                !--- Convert sea ice thickness (meters) to cell average sea ice water
                !--- equivalent (kg m-2) by multiply by sea ice density (kg m-3) and
                !--- sea ice fraction.
                !--- The constant sea ice density used in NEMO/LIM2 is rhoic = 900 kg m-3
                !--- OIceFrc will have been transfered from the ocean along with OIceTck
                !--- so both these fields will be present for the current time (store_time)
                !xxx nwrds = store_readRecord("OIceFrc", ocn_gid, wrk2, time=store_time)
                nwrds = store_readRecord("OIceFrc", ocn_gid, wrk2, idx=-1)
                call  ReadRecord_agcm_err("OIceFrc", -1)
                sea_ice_density = 900.0
                wrk1 = wrk1 * wrk2 * sea_ice_density

              case ("weighted ice and snow")
                !--- OIceTck is the cell averaged sea ice thickness (fr_i*ht_i in NEMO)
                !--- Convert sea ice thickness (meters) to sea ice water equivalent (kg m-2)
                !--- by multiply by sea ice density (kg m-3)
                !--- The constant sea ice density used in NEMO/LIM2 is rhoic = 900 kg m-3
                sea_ice_density = 900.0
                wrk1 = wrk1 * sea_ice_density

              case ("weighted iwe and swe")
                !--- OIceTck is the sea ice water equivalent (rhoic*fr_i*ht_i in NEMO)
                !--- This is what SIC_atm should be. No processing is required.

              case default
                write(6,*)"prepare_agcm_fields: Invalid description for OIceTck --> ", &
                          trim(nemo_namsbc_cpl_cldes(3)),' <--'
                call flush(6)
                call xit("prepare_agcm_fields",-3)

            end select

            src_fld_name = "OIceTck-IWE"

#ifdef DBG_bfb
call dump_farray("OIceTck_before_regrid",wrk1,idx=cpl_elapsed_time_secs)
#endif

            !--- Interpolate this onto the ATM grid
            call local_regrid_Field("OIceTck", wrk1, "SIC_atm", wrk2)

#ifdef DBG_bfb
call dump_farray("SIC_atm_after_regrid",wrk2,idx=cpl_elapsed_time_secs)
#endif

            !--- Interpolation can lead to small negative values
            where ( wrk2 < 0.0 ) wrk2 = 0.0

          case ("SNO_atm")
            !--- Interpolated from NEMO variable OSnwTck
            !--- if sn_snd_thick%cldes = 'ice and snow' then
            !---     OSnwTck = Snow thickness
            !--- if sn_snd_thick%cldes = 'weighted ice and snow' then
            !---     OSnwTck = Cell average snow thickness
            !--- if sn_snd_thick%cldes = 'weighted iwe and swe'
            !---     OSnwTck = Snow water equivalent

            !--- Get field data from internal storage
            !xxx nwrds = store_readRecord("OSnwTck", ocn_gid, wrk1, time=store_time)
            nwrds = store_readRecord("OSnwTck", ocn_gid, wrk1, idx=-1)
            call  ReadRecord_agcm_err("OSnwTck", -1)
write(6,*)"SNO_atm: Received OSnwTck from NEMO  min=",minval(wrk1),"  max=",maxval(wrk1)
write(6,*)"SNO_atm: desc = ",trim(adjustl(nemo_namsbc_cpl_cldes(3)))

            !--- nemo_namsbc_cpl_cldes(3) is sn_snd_thick%cldes sent from NEMO
            select case ( trim(adjustl(nemo_namsbc_cpl_cldes(3))) )
              case ("ice and snow")
                !--- OSnwTck is the non-averaged snow thickness over sea ice (hs_i in NEMO)
                !--- Convert snow thickness (meters) to cell average snow water
                !--- equivalent (kg m-2) by multiply by snow density (kg m-3) and
                !--- sea ice fraction.
                !--- The constant snow density used in NEMO/LIM2 is rhosn = 330 kg m-3
                !--- OIceFrc will have been transfered from the ocean along with OSnwTck
                !--- so both these fields will be present for the current time (store_time)
                !xxx nwrds = store_readRecord("OIceFrc", ocn_gid, wrk2, time=store_time)
                nwrds = store_readRecord("OIceFrc", ocn_gid, wrk2, idx=-1)
                call  ReadRecord_agcm_err("OIceFrc", -1)
                snow_density = 330.0
                wrk1 = wrk1 * wrk2 * snow_density

                !--- Interpolate this onto the ATM grid
                call local_regrid_Field("OSnwTck", wrk1, "SNO_atm", wrk2)

              case ("weighted ice and snow")
                !--- OSnwTck is the cell averaged snow thickness over sea ice (fr_i*hs_i in NEMO)
                !--- Convert snow thickness (meters) to snow water equivalent (kg m-2)
                !--- by multiply by snow density (kg m-3)
                !--- The constant snow density used in NEMO/LIM2 is rhosn = 330 kg m-3
                snow_density = 330.0
                wrk1 = wrk1 * snow_density

                !--- Interpolate this onto the ATM grid
                call local_regrid_Field("OSnwTck", wrk1, "SNO_atm", wrk2)

              case ("weighted iwe and swe")
                !--- OSnwTck is the snow water equivalent over sea ice (rhosn*fr_i*hs_i in NEMO)
                !--- The AGCM wants snow mass per unit area over sea ice fraction (ie not cell average)
                !--- So we need to divide by OIceFrc where it is non-zero
                nwrds = store_readRecord("OIceFrc", ocn_gid, wrk3, idx=-1)
                call  ReadRecord_agcm_err("OIceFrc", -1)

                !--- Interpolate this onto the ATM grid, interpolating numerator
                !--- and denominator prior to dividing by ice fraction
                call local_regrid_Field("OSnwTck", wrk1, "SNO_atm", wrk2)
                call local_regrid_Field("OIceFrc", wrk3, "SNO_atm", wrk4)
                where (wrk4 >  0.0_8) wrk2 = wrk2/wrk4
                where (wrk4 <= 0.0_8) wrk2 = 0.0_8

              case default
                write(6,*)"prepare_agcm_fields: Invalid description for OSnwTck --> ", &
                          trim(nemo_namsbc_cpl_cldes(3)),' <--'
                call flush(6)
                call xit("prepare_agcm_fields",-4)

            end select

            src_fld_name = "OSnwTck-SWE"

            !--- Interpolation can lead to small negative values
            where ( wrk2 < 0.0 ) wrk2 = 0.0
write(6,*)"SNO_atm: Sending SNO_atm to AGCM  min=",minval(wrk2),"  max=",maxval(wrk2)

          case ("CO2flx_atm")
            !--- CO2 flux received from NEMO and sent to the AGCM
            !--- Interpolated from NEMO variable O_CO2FLX
            nwrds = store_readRecord("O_CO2FLX", ocn_gid, wrk1, idx=-1)
            call  ReadRecord_agcm_err("O_CO2FLX", -1)

write(6,*)"O_CO2FLX: Received from NEMO  min=",minval(wrk1),"  max=",maxval(wrk1)

            !--- Convert from (mol C m-2 s-1) to (kg CO2 m-2 s-1), and switch to +ve up (*-1)
            wrk1 = wrk1*0.044011*-1.0

            !--- Interpolate this onto the ATM grid
            call local_regrid_Field("O_CO2FLX", wrk1, "CO2flx_atm", wrk2)

write(6,*)"CO2flx_atm: Sent to AGCM  min=",minval(wrk2),"  max=",maxval(wrk2)
call flush(6)

          case ("RES_atm")
            !--- Residual over lakes
            !--- This is not used when coupled with NEMO

            !--- Send a zero field
            if ( associated(wrk2) ) deallocate(wrk2)
            allocate( wrk2(sz_atm_data) )
            wrk2(:) = 0.0

          case default
            !--- Only the field names listed above are currently recognized
            !--- Add additional cases as required
            write(6,'(3a)') &
              'prepare_agcm_fields: Unknown field name --> ',trim(cpl2atm_name(idx1)),' <--'
            call flush(6)
            call xit("prepare_agcm_fields",-90)
        end select

        if ( verbose > 5 ) then
          call print_var_stats(wrk2, name=cpl2atm_name(idx1), &
                               pfx="prepare_agcm_fields: "//trim(cpl2atm_name(idx1)))
        endif

        !--- Append this data to internal storage
        !--- These fields are associated with the ocn group because they were
        !--- interpolated from data that was received from the ocean
        gid = filter_group_ID(cpl2atm_name(idx1), ocn_gid)
        call store_AppendRecord(cpl2atm_name(idx1), gid, wrk2, create=.true.)

      enddo

      !--- Clean up
      if ( associated(wrk1) ) deallocate(wrk1)
      if ( associated(wrk2) ) deallocate(wrk2)
      if ( associated(wrk3) ) deallocate(wrk3)
      if ( associated(wrk4) ) deallocate(wrk4)

      contains

        subroutine local_regrid_field(src_name, src_data, dst_name, dst_data)
          character(*) :: src_name, dst_name
          real(kind=8), intent(inout), pointer :: src_data(:), dst_data(:)

          if ( len_trim(src_fld_name) < 1 ) then
            src_fld_name = trim(src_name)
          endif

          !--- Conservative-fracarea remap from NEMO to AGCM grid with land values on
          !--- the NEMO grid masked (land does not take part in the remap calculation)
          call regrid_Field(src_name, src_data, dst_name, dst_data, &
                            method="conservative", norm_type="fracarea", &
                            src_fld_name=trim(src_fld_name), &
                            src_mask_values=(/0/), missing_value=0.0_8)

        end subroutine local_regrid_field

        subroutine ReadRecord_agcm_err(name, err)
          character(*) :: name
          integer :: err
          if ( nwrds < 0 ) then
            !--- The field was not found in storage
            write(6,*)"prepare_agcm_fields: ",trim(cpl2atm_name(idx1)), &
                " ...",trim(name)," is missing from storage."
            call xit("prepare_agcm_fields",err)
          endif
        end subroutine ReadRecord_agcm_err

    end subroutine prepare_agcm_fields

    subroutine replace_field(name, nt, wrk)
      character(*), intent(in)    :: name
      integer,      intent(in)    :: nt
      real(kind=8), intent(inout) :: wrk(:)

      !--- Local
      integer(kind=4) :: fid_hist = -1
      integer :: nx_hist, ny_hist
      real(kind=8), pointer :: hist_data(:,:)
      character(256) :: hist_file

      nullify(hist_data)

      write(6,*)"replace_field: ",trim(name),"   nt=",nt," cpl_elapsed_time_secs=", &
                  cpl_elapsed_time_secs
      call flush(6)

      hist_file = " "
      hist_file = "/fs/dev/crb/had02/data/bcpl21__t462_2451_m01_saved_hist.001"

      !--- Open the history file
      call ncdf_open(fid_hist, trim(hist_file), mode="read")

      !--- Read data for the variable at time slice nt
      call ncdf_read_var(fid_hist, trim(name), hist_data, nt=nt)

      !--- Close the history file
      call ncdf_close(fid_hist)

      !--- Verify dimension sizes
      nx_hist = size(hist_data,dim=1)
      ny_hist = size(hist_data,dim=2)
      if ( nx_hist*ny_hist /= size(wrk) ) then
        write(6,*)"replace_field: Size miss match reading ",trim(name)," from ",trim(hist_file)
        write(6,*)"               nx_hist, ny_hist, size(wrk): ",nx_hist,ny_hist,size(wrk)
        call flush(6)
        call xit("replace_field",-1)
      endif

      write(6,*)"replace_field: input wrk array min, max, avg: ", &
                minval(wrk),maxval(wrk),sum(wrk)/real(size(wrk),kind=8)
      write(6,*)"replace_field: local hist_data min, max, avg: ", &
                minval(hist_data),maxval(hist_data),sum(hist_data)/real(size(hist_data),kind=8)

      !--- Overwrite the user suppled wrk array
      wrk = reshape(hist_data, (/nx_hist*ny_hist/))
      deallocate( hist_data )

      write(6,*)"replace_field: output ",trim(name),"   min, max, avg: ", &
                minval(wrk),maxval(wrk),sum(wrk)/real(size(wrk),kind=8)
      call flush(6)

    end subroutine replace_field

end module field_ops
