module history
  !--- History file related procedures for CanCPL

  use com_cpl
  use ccc_esmf
  use ccc_netcdf
  use store

  implicit none

  private
  public :: history_write

  !--- Define grid sizes used in the history file using info from com_cpl
  !--- nlon_a and nlon_o DO NOT include overlap longitude
  integer(kind=4), parameter :: nlon_atm = nlon_a+olap_a  !--- include overlap lon
  integer(kind=4), parameter :: nlat_atm = nlat_a

  integer(kind=4), parameter :: nlon_ocn = nlon_o+olap_o  !--- include overlap lon
  integer(kind=4), parameter :: nlat_ocn = nlat_o

  !--- A derived type to hold meta info about variables written to the history file
  type hist_var_t
    !--- name is a "known" name found in the current cpl_var_list
    character(len=32) :: name = " "
    !--- Group ID that is associated with this variable
    integer :: gid = 0
    !--- The name of the grid this variable is defined on (e.g. atm, ocn, ...)
    character(128) :: grid
    !--- Dimension sizes of the 2D variable to be written to the history file
    integer ::  nx = 0
    integer ::  ny = 0
    !--- Indicate type of variable (either prognostic or flux)
    logical :: isprog = .false.
    logical :: isflux = .false.
    !--- Pointer to the next list element
    type(hist_var_t), pointer :: next => null()
  end type hist_var_t

  type(hist_var_t), save, pointer :: hist_var_list => null()

  logical, save :: history_initialized = .false.

  integer(kind=8), parameter :: secs_per_day=86400
  integer(kind=8), parameter :: secs_per_year=31536000

  !--- A list of variable names for variables that are defined as type prognostic
  character(len=32), dimension(5), parameter ::  prog_names = &
    (/ "SICN_atm                        ", &
       "SIC_atm                         ", &
       "SNO_atm                         ", &
       "GT_atm                          ", &
       "CO2flx_atm                      " /)

  contains

!xxx    function is_var_prog(name) result(isprog)
!xxx      character(*), intent(in) :: name
!xxx      !--- Local
!xxx      logical :: isprog
!xxx      integer :: i
!xxx      isprog = .false.
!xxx      do i=1,len(prog_names)
!xxx        if (trim(adjustl(prog_names(i))) .eq. trim(adjustl(name))) then
!xxx          isprog = .true.
!xxx          exit
!xxx        endif
!xxx      enddo
!xxx      if (isprog) then
!xxx        write(6,*)"is_var_prog: name=",trim(name),"  is a prognostic variable"
!xxx      else
!xxx        write(6,*)"is_var_prog: name=",trim(name),"  is a flux variable"
!xxx      endif
!xxx      call flush(6)
!xxx    end function is_var_prog

    !***************************************************************************
    !--- Determine if a variable originates from the ocean (isprog=T) or from
    !--- the atm (isprog=F)
    !***************************************************************************
    function is_var_prog(name) result(isprog)
      character(*), intent(in) :: name
      !--- Local
      logical :: isprog

      integer :: idx1, idx2, vtype
      character(32) :: get_model, put_model, get_var, put_var
      integer :: verbose=1
      logical :: found

      if ( .not. associated(event) ) then
        write(6,*)"is_var_prog: event is not defined."
        write(6,*)"              Unable to determine variable origin."
        call xit("is_var_prog",-1)
      endif

      if ( nevents < 1 ) then
         write(6,*)"is_var_prog: nevents is out of range. nevents=",nevents
         call xit("is_var_prog",-2)
      endif

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

        found=.false.
        vtype=0
        get_model=" "
        get_var  =" "
        put_model=" "
        put_var  =" "
        get_model=trim(adjustl(event(idx1)%get_model))
        get_var  =trim(adjustl(event(idx1)%get_var))
        put_model=trim(adjustl(event(idx1)%put_model))
        put_var  =trim(adjustl(event(idx1)%put_var))

        !--- Variables to be sent to the ocean
        if ( trim(put_model) .eq. "ocn" ) then
          if (trim(put_var) .eq. trim(adjustl(name))) then
            !--- The input variable (name) is to be sent to the ocean
            found = .true.
            !--- These are defined as type flux, meaning they are derived
            !--- from agcm (flux) variables
            vtype=1
          endif
        endif

        !--- Variables to be sent to the agcm
        if ( trim(put_model) .eq. "atm" ) then
          if (trim(put_var) .eq. trim(adjustl(name))) then
            !--- The input variable (name) is to be sent to the agcm
            found = .true.
            !--- These are defined as type prognostic, meaning they are derived
            !--- from ocean variables
            vtype=2
          endif
        endif

        !--- Variables to be sent from the ocean
        if ( trim(get_model) .eq. "ocn" ) then
          if (trim(get_var) .eq. trim(adjustl(name))) then
            !--- The input variable (name) is to be sent from the ocean
            found = .true.
            vtype=3
          endif
        endif

        !--- Variables to be sent from the agcm
        if ( trim(get_model) .eq. "atm" ) then
          if (trim(get_var) .eq. trim(adjustl(name))) then
            !--- The input variable (name) is to be sent from the agcm
            found = .true.
            vtype=4
          endif
        endif

        if ( verbose > 2 ) then
          write(6,'(3a,i3,2a)') "is_var_prog: name=",trim(name)," event=",idx1, &
            " action=",trim(event(idx1)%action)
          write(6,'(3a,i3,4a)') "is_var_prog: name=",trim(name)," event=",idx1, &
            " get_var=",trim(get_var)," get_model=",trim(get_model)
          write(6,'(3a,i3,4a)') "is_var_prog: name=",trim(name)," event=",idx1, &
            " put_var=",trim(put_var)," put_model=",trim(put_model)
          write(6,'(3a,i3,a,i2,a,l1)') "is_var_prog: name=",trim(name)," event=",idx1, &
            " vtype=",vtype," found=",found
          call flush(6)
        endif

        !--- Accept the first hit
        if ( found ) exit

      enddo

      !--- Assign isprog according to vtype
      isprog=.false.
      if ( vtype == 2 ) isprog=.true.

      if ( verbose > 0 ) then
        write(6,'(3a,i2,a,l1)') &
          "========== is_var_prog: var name=",trim(name)," vtype=",vtype," isprog=",isprog
        call flush(6)
      endif

    end function is_var_prog

    !***************************************************************************
    !--- Copy one hist_var_t variable to another
    !***************************************************************************
    subroutine copy_hist_var_t(var_in, var_out)
      type(hist_var_t), intent(in) :: var_in
      type(hist_var_t), pointer, intent(inout) :: var_out
      nullify(var_out)
      allocate(var_out)
      var_out%name   = var_in%name
      var_out%gid    = var_in%gid
      var_out%grid   = var_in%grid
      var_out%nx     = var_in%nx
      var_out%ny     = var_in%ny
      var_out%isprog = var_in%isprog
      var_out%isflux = var_in%isflux
      if ( associated(var_in%next) ) then
        var_out%next => var_in%next
      else
        var_out%next => null()
      endif
    end subroutine copy_hist_var_t

    !***************************************************************************
    !--- Print info about a single hist_var_t variable
    !***************************************************************************
    subroutine print_hist_var_t(var)
      type(hist_var_t) :: var
      write(6,*)"  name=",trim(var%name)
      write(6,*)"   gid=",var%gid
      write(6,*)"  grid=",trim(var%grid)
      write(6,*)"    nx=",var%nx
      write(6,*)"    ny=",var%ny
      write(6,*)"isprog=",var%isprog
      write(6,*)"isflux=",var%isflux
      call flush(6)
    end subroutine print_hist_var_t

    !***************************************************************************
    !--- Append a history variable to the internal list
    !***************************************************************************
    subroutine hist_var_list_append(new_var)
      type(hist_var_t), intent(in) :: new_var
      type(hist_var_t), pointer :: curr_var
      if ( .not. associated(hist_var_list) ) then
        !--- The list is empty, assign the first element
        call copy_hist_var_t(new_var, hist_var_list)
        hist_var_list%next => null()
        return
      endif
      !--- Find the last element in the list then append
      nullify(curr_var)
      curr_var => hist_var_list
      do
        if ( associated(curr_var%next) ) then
          curr_var => curr_var%next
          cycle
        else
          !--- End of list
          !--- Append the new variable and exit loop
          call copy_hist_var_t(new_var, curr_var%next)
          curr_var => curr_var%next
          curr_var%next => null()
          exit
        endif
      enddo
    end subroutine hist_var_list_append

    !***************************************************************************
    !--- Print info about all history variables in the internal list
    !***************************************************************************
    subroutine hist_var_list_print()
      type(hist_var_t), pointer :: curr_var
      integer :: nvar
      if ( .not. associated(hist_var_list) ) then
        !--- The list is empty
        write(6,*)"hist_var_list_print: The history variable list is empty."
        call flush(6)
        return
      endif
      nullify(curr_var)
      curr_var => hist_var_list
      nvar = 0
      do
        if ( associated(curr_var) ) then
          nvar = nvar+1
          write(6,*)"=== history variable ",nvar
          call print_hist_var_t(curr_var)
        endif
        if ( associated(curr_var%next) ) then
          curr_var => curr_var%next
          cycle
        endif
        exit
      enddo
    end subroutine hist_var_list_print

    !***************************************************************************
    !--- Determine if a history variable currently exists in the internal list
    !***************************************************************************
    function hist_var_list_find(name, gid) result(found)
      character(*), intent(in) :: name
      integer,      intent(in) :: gid
      logical :: found

      !--- Local
      type(hist_var_t), pointer :: curr_var

      if ( .not. associated(hist_var_list) ) then
        !--- The list is empty
        found = .false.
        return
      endif
      found = .false.
      nullify(curr_var)
      curr_var => hist_var_list
      do
        if ( associated(curr_var) ) then
          if ( (trim(adjustl(name)) .eq. &
                trim(adjustl(curr_var%name))) .and. &
               (gid == curr_var%gid) ) then
            found = .true.
            exit
          endif
        endif
        if ( associated(curr_var%next) ) then
          curr_var => curr_var%next
          cycle
        endif
        exit
      enddo
    end function hist_var_list_find

    !***************************************************************************
    !--- Add a history variable to the internal list
    !***************************************************************************
    subroutine register_hist_var(name, gid, grid)
      !--- The name of the variable
      character(*), intent(in) :: name

      !--- The group ID for this variable
      integer,      intent(in) :: gid

      !--- A string to identify the grid this variable is defined on
      !--- This is only required when the input variable name is not defined
      !--- in the global variable cpl_var_list
      character(*), intent(in), optional :: grid

      !--- Local
      type(hist_var_t)  :: curr_var
      type(cpl_vinfo_t), pointer :: curr_vinfo
      integer :: status

      if ( hist_var_list_find(name, gid) ) then
        !--- This variable already exists in the list
        write(6,*)"register_hist_var: Variable ",trim(name),"  gid=",gid," exits."
        call flush(6)
        return
      endif

      nullify(curr_vinfo)
      allocate(curr_vinfo)
      curr_vinfo = find_cpl_vinfo( name=trim(adjustl(name)), status=status )
      if ( status /= 0 ) then
        !--- The record does not exist in cpl_var_list
        !--- If grid is present then use it to create the cpl_var_list entry
        !--- If grid is not present then issue a warning and return (do not abort)
        if ( present(grid) ) then
          !--- Create an entry in cpl_var_list for this variable
          curr_vinfo%name     = trim(adjustl(name))
          curr_vinfo%grid     = trim(adjustl(grid))
          curr_vinfo%tag      = cpl_var_new_tag(6000)
          curr_vinfo%ncopy    = 1
          curr_vinfo%ccc_name = "XXX"
          curr_vinfo%regrid%method   = "conservative"
          curr_vinfo%regrid%norm_opt = "fracarea"
          curr_vinfo%regrid%masked   = .false.
          select case( trim(adjustl(grid)) )
            case( "atm" )
              curr_vinfo%freq     = cpl_atm_freq
              curr_vinfo%size     = (nlon_a+olap_a)*nlat_a
              curr_vinfo%nlon     = nlon_a
              curr_vinfo%nlat     = nlat_a
              curr_vinfo%olap     = olap_a
            case( "ocn" )
              curr_vinfo%freq     = cpl_ocn_freq
              curr_vinfo%size     = (nlon_o+olap_o)*nlat_o
              curr_vinfo%nlon     = nlon_o
              curr_vinfo%nlat     = nlat_o
              curr_vinfo%olap     = olap_o
            case( "msg" )
              curr_vinfo%freq     = -1
              curr_vinfo%size     = 8
              curr_vinfo%nlon     = 8
              curr_vinfo%nlat     = 1
              curr_vinfo%olap     = 0
            case default
              write(6,*)"register_hist_var: Invalid grid --> ",trim(grid)," <--"
              call xit("register_hist_var",-1)
          end select
          !--- Add this variable defintion to cpl_var_list
          call add_cpl_var(curr_vinfo, exists=0)
        else
          !--- Abort when the entry cannot be found or created
          write(6,'(3a)')"register_hist_var: No info for ",trim(name), &
                         " was found. It will not be written to the history file."
          call flush(6)
          return
        endif
      endif

      !--- Define this history variable and add it to the list
      curr_var%name = trim(name)
      curr_var%gid  = gid
      curr_var%grid = trim(adjustl(curr_vinfo%grid))
      curr_var%nx   = curr_vinfo%nlon + curr_vinfo%olap
      curr_var%ny   = curr_vinfo%nlat
      curr_var%isprog = is_var_prog(name)
      curr_var%isflux = .not. curr_var%isprog
      curr_var%next => null()
      call hist_var_list_append(curr_var)

      deallocate(curr_vinfo)

    end subroutine register_hist_var

    !***************************************************************************
    !--- Initialize the list of history variables (aka hist_var_list)
    !***************************************************************************
    subroutine history_init()
      integer :: idx1, idx2, gid
      character(32) :: name
      integer :: verbose=3
      logical :: found
      type(hist_var_t) :: new_var

      !--- Do nothing if history_init has already been called
      if ( history_initialized ) return

      if ( .not. associated(event) ) then
        write(6,*)"history_init: event is not defined."
        write(6,*)"              Unable to initialize history variable list."
        call xit("history_init",-1)
      endif

      if ( nevents < 1 ) then
         write(6,*)"history_init: nevents is out of range. nevents=",nevents
         call xit("history_init",-2)
      endif

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

        !--- Variables to be sent to the ocean
        if ( trim(adjustl(event(idx1)%put_model)) .eq. "ocn" ) then
          name = trim(adjustl(event(idx1)%put_var))
          gid  = filter_group_ID(event(idx1)%put_var, atm_gid)
          call register_hist_var(name, gid)
        endif

        !--- Variables to be sent from the ocean
        if ( trim(adjustl(event(idx1)%get_model)) .eq. "ocn" ) then
          name = trim(adjustl(event(idx1)%get_var))
          gid  = filter_group_ID(event(idx1)%get_var, ocn_gid)
          call register_hist_var(name, gid)
        endif

        !--- Variables to be sent to the agcm
        if ( trim(adjustl(event(idx1)%put_model)) .eq. "atm" ) then
          name = trim(adjustl(event(idx1)%put_var))
          if ( atm_forcing_from_file ) then
            !--- Origin is atm when reading AGCM forcing data from a file
            gid  = filter_group_ID(event(idx1)%put_var, atm_gid)
          else
            gid  = filter_group_ID(event(idx1)%put_var, ocn_gid)
          endif
          call register_hist_var(name, gid)
        endif

        !--- Variables to be sent from the agcm
        if ( trim(adjustl(event(idx1)%get_model)) .eq. "atm" ) then
          name = trim(adjustl(event(idx1)%get_var))
          gid  = filter_group_ID(event(idx1)%get_var, atm_gid)
          call register_hist_var(name, gid)
        endif

      enddo

      !--- Add variables that are used only in the coupler

#ifdef bulk_in_cpl
      call register_hist_var("O_dQlwdT", atm_gid)
      call register_hist_var("O_dQshdT", atm_gid)
      call register_hist_var("O_dQlhdT", atm_gid)
#endif
      if ( found_in_event_list("GTICE_atm") ) then
        !--- GTICE_atm will be interpolated onto GTICE_ocn then written to the history file
        call register_hist_var("GTICE_ocn", atm_gid, grid='ocn')
      else
        write(6,'(2a)') &
          "history_init: GTICE_atm is not in the event list.", &
          " GTICE_ocn will not be written to the history file."
        call flush(6)
      endif
      if ( found_in_event_list("O_TepIce") ) then
        !--- O_TepIce will be interpolated onto O_TepIce_atm then written to the history file
        call register_hist_var("O_TepIce_atm", ocn_gid, grid='atm')
      else
        write(6,'(2a)') &
          "history_init: O_TepIce is not in the event list.", &
          " O_TepIce_atm will not be written to the history file."
        call flush(6)
      endif

      if ( verbose > 2 ) then
        write(6,*)"Printing coupler history variable list"
        call flush(6)
        call hist_var_list_print()
      endif

      history_initialized = .true.

    end subroutine history_init

    !***************************************************************************
    !--- Determine if a user supplied variable name exists in the event list
    !***************************************************************************
    function found_in_event_list(vname) result(found)
      !--- The name of a variable that could be in the event list
      character(*), intent(in) :: vname

      !--- Local
      logical :: found
      integer :: idx1
      integer :: verbose = 1

      if ( .not. associated(event) ) then
        write(6,*)"found_in_event_list: The event list is not defined."
        write(6,*)"              Unable to determine if ",trim(vname)," exists."
        call xit("found_in_event_list",-1)
      endif

      if ( nevents < 1 ) then
         write(6,*)"found_in_event_list: nevents is out of range. nevents=",nevents
         call xit("found_in_event_list",-2)
      endif

      found=.false.
      do idx1=1,nevents
        if ( event(idx1)%is_msg ) then
          !--- Ignore message variables
          cycle
        endif

        if ( verbose > 2 ) then
          write(6,'(a,i5,6a)') &
            "found_in_event_list: idx=",idx1,"  vname=",trim(vname), &
            "  get_var=",trim(event(idx1)%get_var), &
            "  put_var=",trim(event(idx1)%put_var)
          call flush(6)
        endif

        if ( trim(adjustl(event(idx1)%get_var)) .eq. trim(adjustl(vname)) ) then
          !--- This variable name is sent from event(idx1)%get_model
          found=.true.
          if ( verbose > 1 ) then
            write(6,'(5a)') &
              "found_in_event_list: Variable name --> ",trim(vname)," <-- ", &
              "is sent by ",trim(event(idx1)%get_model)
            call flush(6)
          endif
        endif

        if ( trim(adjustl(event(idx1)%put_var)) .eq. trim(adjustl(vname)) ) then
          !--- This variable name is received by event(idx1)%put_model
          found=.true.
          if ( verbose > 1 ) then
            write(6,'(5a)') &
              "found_in_event_list: Variable name --> ",trim(vname)," <-- ", &
              "is received by ",trim(event(idx1)%put_model)
            call flush(6)
          endif
        endif

        if (found) then
          if ( verbose > 0 ) then
            write(6,'(3a)') &
              "found_in_event_list: Found variable --> ",trim(vname)," <-- "
            call flush(6)
          endif
          !--- Accept the first hit
          exit
        endif

      enddo

      if ( .not.found ) then
        if ( verbose > 0 ) then
          write(6,'(3a)') &
            "*WW* found_in_event_list: Did not find variable --> ",trim(vname)," <-- "
          call flush(6)
        endif
      endif

    end function found_in_event_list

    !***************************************************************************
    !--- Write a CanCPL history file
    !***************************************************************************
    subroutine history_write(hist_file, cpl_clock, write_prog, write_flux)
      !--- The name of the history file to be written
      character(*), intent(in) :: hist_file

      !--- The coupler clock
      type(ESMF_Clock), intent(in) :: cpl_clock

      logical, intent(in), optional :: write_prog, write_flux

      !--- Local

      !--- Coupler start time as an ESMF time object
      type(ESMF_Time) :: cpl_start_time
      !--- Coupler stop time as an ESMF time object
      type(ESMF_Time) :: cpl_stop_time
      !--- Coupler current time as an ESMF time object
      type(ESMF_Time) :: cpl_curr_time
      !--- Next time on coupler clock as an ESMF time object
      type(ESMF_Time) :: cpl_next_time
      !--- Reference time on coupler clock as an ESMF time object
      type(ESMF_Time) :: cpl_ref_time
      !--- Simulation time on coupler clock as an ESMF time object
      type(ESMF_TimeInterval) :: cpl_sim_time
      !--- The coupler time step
      type(ESMF_TimeInterval) :: cpl_time_step
      !--- The number of times the coupler clock has been advanced
      integer(ESMF_KIND_I8) :: cpl_advanceCount
      !--- ESMF return status
      integer(ESMF_KIND_I) :: rc_esmf
      type(cpl_vinfo_t) :: curr_vinfo
      type(hist_var_t), pointer :: curr_var
      character(len=32) :: time_string
      integer(kind=4) :: fid, ncstat, vid, time_id, rec_dimid
      integer(kind=4) :: start_idx(3), count_idx(3), ntime, time_idx
      integer :: idx1, nwrds, status, gid
      character(32)  :: var_name
      character(64)  :: msg
      character(128) :: cs_file
      character(len=128) :: grid_type
      !--- Space for data that is read from storage and written to the file
      real(kind=8) :: time_val, ptime_val, ftime_val, tmpr8
      integer(kind=8) :: tmpi8
      real(kind=8), pointer :: vals(:) => null()
      real(kind=8), pointer :: vals2d(:,:) => null()
      integer :: verbose = 1
      logical :: isprog, isflux
      !--- The coupler time as a string
      character(20) :: curr_time_string = " "
      character(20) :: next_time_string = " "
      character(20) :: start_time_string = " "
      character(20) :: stop_time_string = " "
      character(20) :: ref_time_string = " "
      character(20) :: sim_time_string = " "
      character(20) :: time_step_string = " "
      real(ESMF_KIND_R8) :: start_elapsed_secs, ref_elapsed_secs
      real(ESMF_KIND_R8) :: curr_elapsed_secs, next_elapsed_secs, stop_elapsed_secs
      real(ESMF_KIND_R8) :: sim_time_secs_r8, time_step_secs_r8
      integer(ESMF_KIND_I) :: month,day_of_month
      integer(ESMF_KIND_I4) :: i4year,i4hour,i4min,i4sec
      integer(ESMF_KIND_I8) :: sim_time_secs_i8, time_step_secs_i8
      integer(kind=impi) :: myrank=0, ierr=0
      integer(kind=8) :: store_time_secs

      !--- Currently the coupler uses only 1 mpi task, but this could change
      call mpi_comm_rank ( MPI_COMM_WORLD, myrank, ierr )

      !--- Get the number of times that the coupler clock has been advanced
      call ESMF_ClockGet(cpl_clock, advanceCount=cpl_advanceCount, rc=rc_esmf)
      call handle_ESMF_err(rc_esmf, "history_write: Problem reading advance count on coupler clock.")
      write(6,'(i3.3,a,i6,a)') &
        myrank,": history_write: The coupler clock has advanced ", cpl_advanceCount, " times."

      !--- Get the start time on the coupler clock
      call ESMF_ClockGet(cpl_clock, startTime=cpl_start_time, rc=rc_esmf)
      call handle_ESMF_err(rc_esmf, "history_write: Problem reading start time on coupler clock.")

      !--- Get the current time on the clock
      call ESMF_ClockGet(cpl_clock, currTime=cpl_curr_time, rc=rc_esmf)
      call handle_ESMF_err(rc_esmf, "history_write: Problem reading current time on coupler clock.")

      !--- Get the stop time on the coupler clock
      call ESMF_ClockGet(cpl_clock, stopTime=cpl_stop_time, rc=rc_esmf)
      call handle_ESMF_err(rc_esmf, "history_write: Problem reading stop time on coupler clock.")

      !--- Get the coupler time step
      call ESMF_ClockGet(cpl_clock, timeStep=cpl_time_step, rc=rc_esmf)
      call handle_ESMF_err(rc_esmf, "history_write: Problem reading timeStep on coupler clock.")

      !--- Get the coupler reference time step
      call ESMF_ClockGet(cpl_clock, refTime=cpl_ref_time, rc=rc_esmf)
      call handle_ESMF_err(rc_esmf, "history_write: Problem reading refTime on coupler clock.")

      !--- Get the coupler simulation time (curr_time - ref_time)
      call ESMF_ClockGet(cpl_clock, currSimTime=cpl_sim_time, rc=rc_esmf)
      call handle_ESMF_err(rc_esmf, "history_write: Problem reading currSimTime on coupler clock.")

      !--- Get what will be the next time on the coupler clock
      call ESMF_ClockGetNextTime(cpl_clock, cpl_next_time, rc=rc_esmf)
      call handle_ESMF_err(rc_esmf, "history_write: Problem reading next time on coupler clock.")

      !--- cpl_start_time was set to 1 coupling interval before the actual start so that the coupler
      !--- clock could be advanced at the top of the time loop in the coupler driver
      !--- Use the coupler reference time in the following comparison rather than coupler start time
      if ( cpl_curr_time == cpl_ref_time ) then
        !--- This is the first time step of this run
        write(6,'(a)')"history_write: The current time is the reference time (start of month)."
        call flush(6)
      endif

      call ESMF_TimeGet(cpl_start_time, yy=i4year, mm=month, dd=day_of_month, h=i4hour, m=i4min, s=i4sec)
      write(6,'(6(a,i5))')"history_write: coupler clock start_time:  year=",i4year,"  month=",month, &
          "  day=",day_of_month,"  hour=",i4hour,"  minute=",i4min,"  sec=",i4sec
      call flush(6)

      call ESMF_TimeGet(cpl_ref_time, yy=i4year, mm=month, dd=day_of_month, h=i4hour, m=i4min, s=i4sec)
      write(6,'(6(a,i5))')"history_write: coupler clock   ref_time:  year=",i4year,"  month=",month, &
          "  day=",day_of_month,"  hour=",i4hour,"  minute=",i4min,"  sec=",i4sec
      call flush(6)

      call ESMF_TimeGet(cpl_curr_time, yy=i4year, mm=month, dd=day_of_month, h=i4hour, m=i4min, s=i4sec)
      write(6,'(6(a,i5))')"history_write: coupler clock  curr_time:  year=",i4year,"  month=",month, &
          "  day=",day_of_month,"  hour=",i4hour,"  minute=",i4min,"  sec=",i4sec
      call flush(6)

      call ESMF_TimeGet(cpl_next_time, yy=i4year, mm=month, dd=day_of_month, h=i4hour, m=i4min, s=i4sec)
      write(6,'(6(a,i5))')"history_write: coupler clock  next_time:  year=",i4year,"  month=",month, &
          "  day=",day_of_month,"  hour=",i4hour,"  minute=",i4min,"  sec=",i4sec
      call flush(6)

      call ESMF_TimeGet(cpl_stop_time, yy=i4year, mm=month, dd=day_of_month, h=i4hour, m=i4min, s=i4sec)
      write(6,'(6(a,i5))')"history_write: coupler clock  stop_time:  year=",i4year,"  month=",month, &
          "  day=",day_of_month,"  hour=",i4hour,"  minute=",i4min,"  sec=",i4sec
      call flush(6)

      !--- Assign the string start_time_string using info read from the coupler clock
      call ESMF_TimeGet(cpl_start_time, s_r8=start_elapsed_secs, timeStringISOFrac=start_time_string)
      do idx1=1,len_trim(start_time_string)
        if ( start_time_string(idx1:idx1) .eq. "T" ) start_time_string(idx1:idx1) = " "
      enddo
      !--- Convert start_elapsed_secs to years
      tmpr8 = start_elapsed_secs/real(secs_per_year,kind=8)
      write(6,'(i3.3,3a,g24.16)') &
        myrank,": history_write: coupler clock start_time=",trim(start_time_string), &
        "  elapsed years: ",tmpr8
      call flush(6)

      !--- Assign the string ref_time_string with the coupler reference time
      call ESMF_TimeGet(cpl_ref_time, s_r8=ref_elapsed_secs, timeStringISOFrac=ref_time_string)
      do idx1=1,len_trim(ref_time_string)
        if ( ref_time_string(idx1:idx1) .eq. "T" ) ref_time_string(idx1:idx1) = " "
      enddo
      !--- Convert ref_elapsed_secs to years
      tmpr8 = ref_elapsed_secs/real(secs_per_year,kind=8)
      write(6,'(i3.3,3a,g24.16)') &
        myrank,": history_write: coupler clock   ref_time=",trim(ref_time_string), &
        "  elapsed years: ",tmpr8
      call flush(6)

      !--- Assign the string curr_time_string with the current coupler time
      call ESMF_TimeGet(cpl_curr_time, s_r8=curr_elapsed_secs, timeStringISOFrac=curr_time_string)
      do idx1=1,len_trim(curr_time_string)
        if ( curr_time_string(idx1:idx1) .eq. "T" ) curr_time_string(idx1:idx1) = " "
      enddo
      !--- Convert curr_elapsed_secs to years
      tmpr8 = curr_elapsed_secs/real(secs_per_year,kind=8)
      write(6,'(i3.3,3a,g24.16)') &
        myrank,": history_write: coupler clock  curr_time=",trim(curr_time_string), &
        "  elapsed years: ",tmpr8
      call flush(6)

      !--- Assign the string stop_time_string with the coupler stop time
      call ESMF_TimeGet(cpl_stop_time, s_r8=stop_elapsed_secs, timeStringISOFrac=stop_time_string)
      do idx1=1,len_trim(stop_time_string)
        if ( stop_time_string(idx1:idx1) .eq. "T" ) stop_time_string(idx1:idx1) = " "
      enddo
      !--- Convert stop_elapsed_secs to years
      tmpr8 = stop_elapsed_secs/real(secs_per_year,kind=8)
      write(6,'(i3.3,3a,g24.16)') &
        myrank,": history_write: coupler clock  stop_time=",trim(stop_time_string), &
        "  elapsed years: ",tmpr8
      call flush(6)

      !--- Assign the string next_time_string with the next coupler time
      call ESMF_TimeGet(cpl_next_time, s_r8=next_elapsed_secs, timeStringISOFrac=next_time_string)
      do idx1=1,len_trim(next_time_string)
        if ( next_time_string(idx1:idx1) .eq. "T" ) next_time_string(idx1:idx1) = " "
      enddo
      !--- Convert next_elapsed_secs to years
      tmpr8 = next_elapsed_secs/real(secs_per_year,kind=8)
      write(6,'(i3.3,3a,g24.16)') &
        myrank,": history_write: coupler clock  next_time=",trim(next_time_string), &
        "  elapsed years: ",tmpr8
      call flush(6)

      !--- Assign the string sim_time_string with the coupler simulation time
      call ESMF_TimeIntervalGet(cpl_sim_time, s_i8=sim_time_secs_i8, s_r8=sim_time_secs_r8, &
                                timeStringISOFrac=sim_time_string)
      do idx1=1,len_trim(sim_time_string)
        if ( sim_time_string(idx1:idx1) .eq. "T" ) sim_time_string(idx1:idx1) = " "
      enddo
      write(6,'(i3.3,3a,g24.16)') &
        myrank,": history_write: sim_time=",trim(sim_time_string),"  elapsed seconds: ",sim_time_secs_r8
      call flush(6)
      write(6,'(a,g24.16)')"history_write: sim_time elapsed seconds: ",sim_time_secs_r8
      write(6,'(a,i10)')   "history_write:    cpl_elapsed_time_secs: ",cpl_elapsed_time_secs
      call flush(6)

      !--- Assign the string sim_time_string with the coupler simulation time
      call ESMF_TimeIntervalGet(cpl_time_step, s_i8=time_step_secs_i8, s_r8=time_step_secs_r8, &
                                timeStringISOFrac=time_step_string)
      do idx1=1,len_trim(time_step_string)
        if ( time_step_string(idx1:idx1) .eq. "T" ) time_step_string(idx1:idx1) = " "
      enddo
      write(6,'(i3.3,3a,g24.16)') &
        myrank,": history_write: time_step=",trim(time_step_string),"  seconds: ",time_step_secs_r8
      call flush(6)
      write(6,'(a,g24.16)')"history_write: time step elapsed seconds: ",time_step_secs_r8
      call flush(6)

!xxx      !--- Find out if the next coupler time step is the last
!xxx      if ( cpl_next_time == cpl_stop_time ) then
!xxx        !--- Do not write anything to the history file
!xxx        write(6,'(i3.3,3a)')myrank,": history_write: Next time is ",trim(next_time_string), &
!xxx                  ". This is the last time step. Nothing is written to the history file."
!xxx        call flush(6)
!xxx        return
!xxx      endif

      !--- Open a history file (or create if missing) and return the netcdf file id
      call history_open(hist_file, fid, cpl_ref_time, cpl_stop_time)

      ncstat = nf90_Inquire(fid, unlimitedDimId=rec_dimid)
      msg = "history_write: retrieving unlimited dim ID"
      call nc_error_handler(fid, ncstat, 1, msg=trim(msg))

      !--- Determine the number of time records currently in the history file
      ncstat = nf90_inq_dimid(fid, "time", time_id)
      call nc_error_handler(fid, ncstat, 1, "retrieving time ID")
      ncstat = nf90_inquire_dimension(fid, time_id, len=ntime)
      call nc_error_handler(fid, ncstat, 2, "retrieving max record length")

      !--- The simulation time in days (current_time - reference_time)
      !--- Note that the the first call will be for 1 coupler time step before the reference time
      !--- and so the sim time determined above will be negative 1 time step on the first call.
      !--- Increment sim time by 1 coupler time step to set the first write time to zero.
      !--- This time is used as the time slice value that is written to the history file.
      tmpi8 = sim_time_secs_i8 + time_step_secs_i8
      time_val = real(tmpi8,kind=8)/real(secs_per_day,kind=8)

      write(6,'(a,i6,a)') &
        "history_write: The coupler clock has advanced ", cpl_advanceCount, " times."
      write(6,'(2a)')"history_write: Current time on coupler clock: ",trim(curr_time_string)
      write(6,'(a,g16.7)')"history_write:       write_time_days=",time_val
      write(6,'(a,i10)')  "history_write:       write_time_secs=",tmpi8
      write(6,'(a,i10)')  "history_write: cpl_elapsed_time_secs=",cpl_elapsed_time_secs
      call flush(6)

      write(6,'(i3.3,a,5i8)') &
        myrank,": history_write: fid, time_id, rec_dimid, ntime, cpl_time_step_iter = ", &
                                 fid, time_id, rec_dimid, ntime, cpl_time_step_iter
      call flush(6)

      if ( ntime > 0 ) then
        time_idx = ntime + 1
      else
        time_idx = 1
      endif
if (.false.) then
      if ( ntime > 0 ) then
        !--- Read the current time values
        if ( associated(vals) ) deallocate(vals)
        allocate( vals(ntime) )
        ncstat = nf90_get_var( fid, time_id, vals, start=(/1_4/), count=(/ntime/) )
        msg = " "
        write(msg,'(a,i4,a,i2)')"history_write: Reading ",ntime," time values from vid=",time_id
        call nc_error_handler(fid, ncstat, 3, msg=trim(msg))

        !--- Set time_idx to the index for the next new time record
        time_idx = ntime + 1

        !--- Look for time_val in the existing time coordinate values
        do idx1=1,size(vals)
          if ( vals(idx1) == time_val ) then
            !--- time_val already exists as a time coordinate
            !--- Reset time_idx to the index of this time value
            time_idx = idx1
            exit
          endif
        enddo
        deallocate(vals)
      else
        time_idx = 1
      endif
endif

     if ( couple_serial ) then
       !--- Serial coupling mode:
       !--- Write the time coordinate value for the next record
       !--- except when the next time step is the last time step
       !--- On this time step flux fields are written but prognostic fields are not
       !--- and since flux fields are always written at the previous time index there will
       !--- be no fields written when cpl_next_time == cpl_stop_time
       if ( cpl_next_time /= cpl_stop_time ) then
         ncstat = nf90_inq_varid(fid, "time", vid)
         call nc_error_handler(fid, ncstat, 4, "getting time variable ID")
         ncstat = nf90_put_var( fid, vid, (/time_val/), start=(/time_idx/), count=(/1_4/) )
         call nc_error_handler(fid, ncstat, 5, "writing time value")
       endif
     else
       !--- Parallel coupling mode:
       !--- Write the time coordinate value for the next record
       ncstat = nf90_inq_varid(fid, "time", vid)
       call nc_error_handler(fid, ncstat, 4, "getting time variable ID")
       ncstat = nf90_put_var( fid, vid, (/time_val/), start=(/time_idx/), count=(/1_4/) )
       call nc_error_handler(fid, ncstat, 5, "writing time value")
     endif

write(6,'(a,g16.8,a)')"history_write: writing at time = ",time_val," days."
call flush(6)

      !--- Define the starting indicies for each prognostic variable write done below
      start_idx(1) = 1
      start_idx(2) = 1
      start_idx(3) = time_idx

      !--- Loop over all entries in the history variable list and write these
      !--- variables to the coupler history file
      nullify(curr_var)
      curr_var => hist_var_list
      do
        if ( associated(curr_var) ) then
          isprog = is_var_prog(curr_var%name)
          isflux = .not.isprog
          if ( cpl_advanceCount == 0 ) then
            !--- This is the first call from the coupler driver (before the main time stepping loop
            !--- has started). The fields written to the history file here will be those read from the
            !--- coupler restart.
            if ( couple_serial ) then
              !--- Serial coupling mode: Write only ocean prognostic variables on the first time step
              if ( isprog ) then
                write(6,'(i3.3,5a)') &
                    myrank,": history_write: First step: coupler clock curr_time=",trim(curr_time_string), &
                    " writing ",trim(curr_var%name)," to the coupler history file."
                call flush(6)
              else
                write(6,'(i3.3,5a)') &
                    myrank,": history_write: First step: coupler clock curr_time=",trim(curr_time_string), &
                    " NOT writing flux variable ",trim(curr_var%name)," to the coupler history file."
                call flush(6)
                curr_var => curr_var%next
                cycle
              endif
            else
              !--- Parallel coupling mode: Write all fields on the first time step
              write(6,'(i3.3,5a)') &
                  myrank,": history_write: First step: coupler clock curr_time=",trim(curr_time_string), &
                  " writing ",trim(curr_var%name)," to the coupler history file."
              call flush(6)
            endif
          endif
          if ( couple_serial ) then
            !--- Serial coupling mode: Do not write ocean prognostic variables on the last time step
            if ( cpl_next_time == cpl_stop_time ) then
              !--- Do not write any prognostic variables to the history file
              if ( isprog ) then
                write(6,'(i3.3,5a)') &
                    myrank,": history_write: Last step: coupler clock curr_time=",trim(curr_time_string), &
                    " NOT writing prognostic variable ",trim(curr_var%name)," to the coupler history file."
                call flush(6)
                curr_var => curr_var%next
                cycle
              endif
            endif
          endif

          !--- Define the number of indicies along each dimension
          count_idx(1) = curr_var%nx
          count_idx(2) = curr_var%ny
          count_idx(3) = 1

          if ( couple_serial ) then
            !--- The time index at which to write depends on the variable type (ie prognostic or flux)
            !--- Serial coupling mode: Offset time by one time step for flux fields
            if ( isprog ) then
              start_idx(3) = time_idx
            else
              start_idx(3) = time_idx - 1
            endif
          else
            !--- Parallel coupling mode: Write all fields at the same time level
            start_idx(3) = time_idx
          endif

          if ( associated(vals2d) ) deallocate(vals2d)
          allocate( vals2d(curr_var%nx,curr_var%ny) )

          !--- Get field data from internal storage
          !--- Get the record with time = cpl_elapsed_time_secs for the current name
!xxx          nwrds = store_readRecord(curr_var%name, curr_var%gid, vals, time=cpl_elapsed_time_secs, &
!xxx                                   time_out=store_time_secs)
          call history_get_field_data(curr_var, vals, nwrds, &
                                      time_in=cpl_elapsed_time_secs, &
                                      time_out=store_time_secs)

          call store_ListRecord(curr_var%name,anygid,pfx="history_write")

          if ( nwrds /= curr_var%nx * curr_var%ny ) then
            write(6,'(3a)')"history_write: Size mismatch reading ", &
                           trim(curr_var%name)," from internal storage."
            write(6,*)"   Expected size ",curr_var%nx*curr_var%ny,"  Found size ",nwrds
            call xit("history_write",-1)
          endif

          !--- Assign these values to a 2d array for the write to the history
          vals2d(:,:) = reshape( vals(1:nwrds), (/curr_var%nx, curr_var%ny/) )

          !--- Write this field to the history file
          ncstat = nf90_inq_varid(fid, trim(curr_var%name), vid)
          if (ncstat /= NF90_NOERR) then
            !--- The current variable does not exists in the history file
            !--- Define it and return the variable ID
            call history_add_var(fid, curr_var%name, curr_var%grid, vid)
          endif
          ncstat = nf90_put_var( fid, vid, vals2d, start=start_idx, count=count_idx )
          call nc_error_handler(fid, ncstat, 2, "writing "//trim(curr_var%name))

          if ( verbose > 0 ) then
            write(6,'(3a,i3,3a,g16.8,a)')"history_write: Wrote ",trim(curr_var%name), &
                           " gid=",curr_var%gid," to ",trim(hist_file),"   time=",time_val," days."
            write(6,'(3a,3i6)')"history_write: ",trim(curr_var%name)," start_idx: ",start_idx(1:3)
            write(6,'(3a,3i6)')"history_write: ",trim(curr_var%name)," count_idx: ",count_idx(1:3)
            write(6,'(3a,i10)')"history_write: ",trim(curr_var%name), &
              " cpl_elapsed_time_secs=",cpl_elapsed_time_secs
            write(6,'(3a,i10)')"history_write: ",trim(curr_var%name), &
              "       store_time_secs=",store_time_secs
            call flush(6)
          endif

          deallocate(vals2d)

          curr_var => curr_var%next
          cycle
        endif
        exit
      enddo
      if ( associated(vals) ) deallocate( vals )

      !--- Close the history file
      ncstat = nf90_close(fid)
      call nc_error_handler(fid, ncstat, 99)

    end subroutine history_write

    !***************************************************************************
    !--- Read data from internal storage (at a particular time)
    !--- given the variable name and variable group ID
    !***************************************************************************
    subroutine history_get_field_data(hist_var, vals, nwrds, time_in, time_out)
      type(hist_var_t), pointer, intent(in) :: hist_var
      real(kind=8), pointer, intent(inout) :: vals(:)
      integer(kind=8), intent(in) :: time_in
      integer(kind=8), intent(out) :: time_out
      integer :: nwrds

      !--- Local
      integer :: local_nwrds
      integer(kind=8) :: local_time_out
      real(kind=8), pointer :: wrk1(:) => null()
      real(kind=8), pointer :: wrk2(:) => null()

      select case(trim(adjustl(hist_var%name)))
      case("GTICE_ocn")
        !--- Interpolate GTICE_ocn from GTICE_atm, if GTICE_atm is available
        if ( store_GetIndex("GTICE_atm", anygid) > 0 ) then
          !--- Read GTICE_atm from storage and interpolate onto GTICE_ocn
          local_nwrds = store_readRecord("GTICE_atm", anygid, wrk1, idx=-1, time_out=time_out)
          if ( .not. cpl_var_name_exists("GTICE_ocn") ) then
             call define_cpl_var_by_name("GTICE_ocn", "ocn")
          endif
          call regrid_atm_to_nemo("GTICE_atm", wrk1, "GTICE_ocn", wrk2)
          nwrds = size(wrk2)
          if (size(wrk2) == hist_var%nx * hist_var%ny) then
            !--- The interpolated field has the expected size
            !--- Allocate and assign the vals array
            if ( associated(vals) ) nullify(vals)
            allocate(vals(hist_var%nx * hist_var%ny))
            vals = wrk2
          else
            write(6,*)"history_get_field_data: Array size mismatch"
            write(6,*)"    Expecting size: ",hist_var%nx * hist_var%ny
            write(6,*)"        Found size: ",size(wrk2)
            call flush(6)
            call xit("history_get_field_data",-1)
          endif
          call store_AppendRecord("GTICE_ocn", atm_gid, vals, create=.true.)
        else
          write(6,*)"history_get_field_data: Requesting GTICE_ocn but GTICE_atm is missing."
          call flush(6)
          call store_ListStorage()
          call xit("history_get_field_data",-2)
        endif
      case("O_TepIce_atm")
        !--- Interpolate O_TepIce_atm from O_TepIce, if O_TepIce is available
        if ( store_GetIndex("O_TepIce", anygid) > 0 ) then
          !--- Read O_TepIce from storage and interpolate onto O_TepIce_atm
          local_nwrds = store_readRecord("O_TepIce", anygid, wrk1, idx=-1, time_out=time_out)
          if ( .not. cpl_var_name_exists("O_TepIce_atm") ) then
             call define_cpl_var_by_name("O_TepIce_atm", "atm")
          endif
          call regrid_nemo_to_atm("O_TepIce", wrk1, "O_TepIce_atm", wrk2)
          nwrds = size(wrk2)
          if (size(wrk2) == hist_var%nx * hist_var%ny) then
            !--- The interpolated field has the expected size
            !--- Allocate and assign the vals array
            if ( associated(vals) ) nullify(vals)
            allocate(vals(hist_var%nx * hist_var%ny))
            vals = wrk2
          else
            write(6,*)"history_get_field_data: Array size mismatch"
            write(6,*)"    Expecting size: ",hist_var%nx * hist_var%ny
            write(6,*)"        Found size: ",size(wrk2)
            call flush(6)
            call xit("history_get_field_data",-3)
          endif
          call store_AppendRecord("O_TepIce_atm", ocn_gid, vals, create=.true.)
        else
          write(6,*)"history_get_field_data: Requesting O_TepIce_atm but O_TepIce is missing."
          call flush(6)
          call store_ListStorage()
          call xit("history_get_field_data",-4)
        endif
      case default
        !--- Get field data from internal storage
        if ( store_GetIndex(hist_var%name, hist_var%gid) > 0 ) then
          !--- Get the record with time = time_in (ie cpl_elapsed_time_secs) for the given name
          nwrds = store_readRecord(hist_var%name, hist_var%gid, vals, time=time_in, time_out=time_out)
        else
          write(6,*)"history_get_field_data: Requested variable ",trim(hist_var%name), &
                    " gid=",hist_var%gid," is missing."
          call flush(6)
          call store_ListStorage()
          call xit("history_get_field_data",-5)
        endif
      end select

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

    end subroutine history_get_field_data

    !***************************************************************************
    !--- Regrid data from the AGCM grid to the NEMO grid
    !--- Mask land values on the destination grid (NEMO grid)
    !***************************************************************************
    subroutine regrid_atm_to_nemo(src_name, src_data, dst_name, dst_data)
      character(*) :: src_name, dst_name
      real(kind=8), intent(inout), pointer :: src_data(:), dst_data(:)

      if ( .not. associated(src_data) ) then
        write(6,'(3a)') &
          "cpl::history::regrid_atm_to_nemo: Unable to interpolate ",trim(src_name), &
          ", missing input data."
        call flush(6)
        call xit("regrid_atm_to_nemo",-1)
      endif

      !--- Conservative-fracarea remap from AGCM to NEMO grid with land values on
      !--- the NEMO grid masked (ie: 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_name), &
                        dst_mask_values=(/0/), missing_value=0.0_8)

    end subroutine regrid_atm_to_nemo

    !***************************************************************************
    !--- Regrid data from the NEMO grid to the AGCM grid
    !--- Mask land values on the source grid (the NEMO grid)
    !***************************************************************************
    subroutine regrid_nemo_to_atm(src_name, src_data, dst_name, dst_data)
      character(*) :: src_name, dst_name
      real(kind=8), intent(inout), pointer :: src_data(:), dst_data(:)

      if ( .not. associated(src_data) ) then
        write(6,'(3a)') &
          "cpl::history::regrid_nemo_to_atm: Unable to interpolate ",trim(src_name), &
          ", missing input data."
        call flush(6)
        call xit("regrid_nemo_to_atm",-1)
      endif

      !--- Conservative-fracarea remap from NEMO to AGCM grid with land values on
      !--- the NEMO grid masked (ie: 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_name), &
                        src_mask_values=(/0/), missing_value=0.0_8)

    end subroutine regrid_nemo_to_atm

    !***************************************************************************
    !--- Add a variable definition to the history file
    !***************************************************************************
    subroutine history_add_var(fid, name, grid, vid)
      integer(kind=4), intent( in) :: fid
      character(*),    intent( in) :: name
      character(*),    intent( in) :: grid
      integer(kind=4), intent(out) :: vid

      !--- Local
      integer (kind=4) :: ncstat, lon_id, lat_id, time_id
      integer :: verbose = 1

      !--- Read dimension IDs
      ncstat = nf90_inq_dimid(fid, "time", time_id)
      call nc_error_handler(fid, ncstat, 1, "history_add_var: retrieving time ID")

      select case( trim(adjustl(grid)) )
        case( "atm" )
          ncstat = nf90_inq_dimid(fid, "lon_a", lon_id)
          call nc_error_handler(fid, ncstat, 1, "history_add_var: retrieving lon ID")
          ncstat = nf90_inq_dimid(fid, "lat_a", lat_id)
          call nc_error_handler(fid, ncstat, 1, "history_add_var: retrieving lat ID")

        case( "ocn" )
          ncstat = nf90_inq_dimid(fid, "nlon_o", lon_id)
          call nc_error_handler(fid, ncstat, 1, "history_add_var: retrieving lon ID")
          ncstat = nf90_inq_dimid(fid, "nlat_o", lat_id)
          call nc_error_handler(fid, ncstat, 1, "history_add_var: retrieving lat ID")

        case default
          write(6,*)"history_add_var: Invalid grid type --> ",trim(grid)," <--"
          call xit("history_add_var",-1)

      end select

      !--- Put the file into define mode
      ncstat = nf90_redef(fid)
      call nc_error_handler(fid, ncstat, 1, &
          "history_add_var: Enter define mode. var="//trim(name))

      !--- Define the variable
      ncstat = nf90_def_var(fid, trim(adjustl(name)), NF90_FLOAT, &
                            (/lon_id, lat_id, time_id/), vid)
      call nc_error_handler(fid, ncstat, 1, "history_add_var: defining "//trim(name))

      !--- Add variable specific meta data
      call nc_add_var_meta_data(name, fid, vid)

      !--- Get out of define mode
      ncstat = nf90_enddef(fid)
      call nc_error_handler(fid, ncstat, 14, &
          "history_add_var: Exit define mode. var="//trim(name))

      if ( verbose > 0 ) then
        write(6,'(5a,4(a,i6))')"history_add_var: Defined var=",trim(name), &
            " on the ",trim(grid)," grid.","  vid=",vid, &
            "  lon_id=",lon_id,"  lat_id=",lat_id,"  time_id=",time_id
        call flush(6)
      endif

    end subroutine history_add_var

    !***************************************************************************
    !--- Open or create a history file
    !***************************************************************************
    subroutine history_open(hist_file, fid, cpl_ref_time, cpl_stop_time)
      character(*),    intent(in)  :: hist_file
      integer(kind=4), intent(out) :: fid

      !--- Coupler reference time as an ESMF time object
      type(ESMF_Time), intent(in) :: cpl_ref_time

      !--- Coupler stop time as an ESMF time object
      type(ESMF_Time), intent(in) :: cpl_stop_time

      !--- Local
      integer (kind=4) :: ncstat, vid, idx
      integer (kind=4) :: time_dim_id, time_id
      integer (kind=4) :: lon_dim_atm_id, lat_dim_atm_id
      integer (kind=4) :: lon_dim_ocn_id, lat_dim_ocn_id
      integer (kind=4) :: lon_atm_id, lat_atm_id
      integer (kind=4) :: lon_ocn_id, lat_ocn_id
      integer (kind=4) :: mask_atm_id, mask_ocn_id
      integer (kind=4) :: atm_fland_id, atm_grid_cell_area_id, ocn_grid_cell_area_id
      integer :: idx1, nvar_atm, nvar_ocn, grid_file_type, verbose=1
      character(32) :: var_name
      logical :: exists, found
      character(64) :: strng
      character(32) :: time_string
      character(8)  :: curr_date
      character(10) :: curr_time

      !--- A temporary variable used to create new entries in cpl_var_list
      type(cpl_vinfo_t), pointer :: tmp_vinfo => null()

      !--- The total number of points in each grid
      integer :: grid_size_atm, grid_size_ocn

      !--- The shape of each grid ie (nx,ny)
      integer :: grid_dims_atm(2), grid_dims_ocn(2)

      !--- lon/lat of grid cell centers
      real(kind=8), pointer ::  mid_lon_atm(:,:), mid_lon_ocn(:,:)
      real(kind=8), pointer ::  mid_lat_atm(:,:), mid_lat_ocn(:,:)

      !--- Integer mask (1 means valid grid cell, 0 means ignore this grid cell)
      integer, pointer :: mask_atm(:,:), mask_ocn(:,:)

      !--- tmp array
      real(kind=8), pointer :: wrk1d(:)
      real(kind=8), pointer :: wrk2d(:,:)
      integer :: nwrds
 
      nullify(tmp_vinfo, mask_atm, mask_ocn, wrk1d, wrk2d)
      nullify(mid_lon_atm, mid_lon_ocn, mid_lat_atm, mid_lat_ocn)

      !--- On return the netcdf file will be open and in data mode

      !--- Ensure the list of history variables is defined
      call history_init()

      inquire( file=trim(hist_file), exist=exists )

      if ( exists ) then
        !--- Open the file for writing
        ncstat = nf90_open(hist_file, NF90_WRITE, fid)
        call nc_error_handler(fid, ncstat, 1, "Opening history file "//trim(hist_file))
        return
      else
        !--- Create the output netcdf file in netcdf 4 format
        ncstat = nf90_create(hist_file, NF90_NETCDF4, fid)
        call nc_error_handler(fid, ncstat, 2, "Creating history file "//trim(hist_file))
      endif

      !--- Write global attributes
      time_string = " "
      call ESMF_TimeGet(cpl_ref_time, timeStringISOFrac=time_string)
      do idx1=1,len_trim(time_string)
        if ( time_string(idx1:idx1) .eq. "T" ) time_string(idx1:idx1) = " "
      enddo
      ncstat = nf90_put_att(fid, NF90_GLOBAL, "model_start", trim(time_string) )
      call nc_error_handler(fid, ncstat, 3, "adding global attribute model_start")

      time_string = " "
      call ESMF_TimeGet(cpl_stop_time, timeStringISOFrac=time_string)
      do idx1=1,len_trim(time_string)
        if ( time_string(idx1:idx1) .eq. "T" ) time_string(idx1:idx1) = " "
      enddo
      ncstat = nf90_put_att(fid, NF90_GLOBAL, "model_stop", trim(time_string) )
      call nc_error_handler(fid, ncstat, 4, "adding global attribute model_stop")

      if ( len_trim(env_runid) > 0 ) then
        ncstat = nf90_put_att(fid, NF90_GLOBAL, "CCCma_runid", trim(env_runid) )
        call nc_error_handler(fid, ncstat, 5, "adding global attribute CCCma_runid")
      endif

      call date_and_time(curr_date, curr_time)
      strng = "created: "//curr_date(1:4)//"-"//curr_date(5:6)//"-"     &
          //curr_date(7:8)//" "//curr_time(1:2)//":"//curr_time(3:4)  &
          //":"//curr_time(5:6)
      ncstat = nf90_put_att(fid, NF90_GLOBAL, "history", trim(strng) )
      call nc_error_handler(fid, ncstat, 6, "adding global attribute history")

      if ( len_trim(env_runid) > 0 ) then
        ncstat = nf90_put_att(fid, NF90_GLOBAL, "CCCma_runid", trim(env_runid) )
        call nc_error_handler(fid, ncstat, 3, "adding global attribute CCCma_runid")
      endif

      if ( len_trim(cpl_build_commit_id) > 0 ) then
        ncstat = nf90_put_att(fid, NF90_GLOBAL, "CanCPL_build_commit", trim(cpl_build_commit_id) )
        call nc_error_handler(fid, ncstat, 7, "adding global attribute CanCPL_build_commit")
      endif

      if ( len_trim(cpl_build_repo_path) > 0 ) then
        ncstat = nf90_put_att(fid, NF90_GLOBAL, "CanCPL_build_repo", trim(cpl_build_repo_path) )
        call nc_error_handler(fid, ncstat, 8, "adding global attribute CanCPL_build_repo")
      endif

      if ( len_trim(cpl_runtime_commit) > 0 ) then
        ncstat = nf90_put_att(fid, NF90_GLOBAL, "CanCPL_runtime_commit", trim(cpl_runtime_commit) )
        call nc_error_handler(fid, ncstat, 7, "adding global attribute CanCPL_runtime_commit")
      endif

      if ( len_trim(cpl_runtime_repo_path) > 0 ) then
        ncstat = nf90_put_att(fid, NF90_GLOBAL, "CanCPL_runtime_repo", trim(cpl_runtime_repo_path) )
        call nc_error_handler(fid, ncstat, 8, "adding global attribute CanCPL_runtime_repo")
      endif

      if ( .not. atm_forcing_from_file ) then
        !--- Do not add any NEMO related meta data when NEMO is not coupled
        if ( len_trim(nemo_commit) > 0 ) then
          ncstat = nf90_put_att(fid, NF90_GLOBAL, "NEMO_build_commit", trim(nemo_commit) )
          call nc_error_handler(fid, ncstat, 7, "adding global attribute NEMO_commit")
        endif

        if ( len_trim(nemo_repo_path) > 0 ) then
          ncstat = nf90_put_att(fid, NF90_GLOBAL, "NEMO_build_repo", trim(nemo_repo_path) )
          call nc_error_handler(fid, ncstat, 8, "adding global attribute NEMO_repo")
        endif

        if ( len_trim(nemo_config) > 0 ) then
          ncstat = nf90_put_att(fid, NF90_GLOBAL, "NEMO_build_config", trim(nemo_config) )
          call nc_error_handler(fid, ncstat, 8, "adding global attribute NEMO_config")
        endif

        if ( len_trim(nemo_build_cancpl_commit) > 0 ) then
          ncstat = nf90_put_att(fid, NF90_GLOBAL, "NEMO_build_cancpl_commit", trim(nemo_build_cancpl_commit) )
          call nc_error_handler(fid, ncstat, 7, "adding global attribute NEMO_build_cancpl_commit")
        endif

        if ( len_trim(nemo_build_cancpl_repo_path) > 0 ) then
          ncstat = nf90_put_att(fid, NF90_GLOBAL, "NEMO_build_cancpl_repo", trim(nemo_build_cancpl_repo_path) )
          call nc_error_handler(fid, ncstat, 8, "adding global attribute NEMO_build_cancpl_repo")
        endif

        if ( len_trim(nemo_runtime_commit) > 0 ) then
          ncstat = nf90_put_att(fid, NF90_GLOBAL, "NEMO_runtime_commit", trim(nemo_runtime_commit) )
          call nc_error_handler(fid, ncstat, 7, "adding global attribute NEMO_runtime_commit")
        endif

        if ( len_trim(nemo_runtime_repo_path) > 0 ) then
          ncstat = nf90_put_att(fid, NF90_GLOBAL, "NEMO_runtime_repo", trim(nemo_runtime_repo_path) )
          call nc_error_handler(fid, ncstat, 8, "adding global attribute NEMO_runtime_repo")
        endif

        if ( len_trim(nemo_runtime_config) > 0 ) then
          ncstat = nf90_put_att(fid, NF90_GLOBAL, "NEMO_runtime_config", trim(nemo_runtime_config) )
          call nc_error_handler(fid, ncstat, 8, "adding global attribute NEMO_runtime_config")
        endif
      endif

      strng = " "
      strng = nf90_inq_libvers()
      do idx1=1,len_trim(strng)
        if ( strng(idx1:idx1) .eq. "$" ) strng(idx1:idx1) = " "
      enddo
      ncstat = nf90_put_att(fid, NF90_GLOBAL, "netcdf_version", trim(strng) )
      call nc_error_handler(fid, ncstat, 9, "adding global attribute netcdf_version")

      !--- Define dimension IDs
      ncstat = nf90_def_dim(fid, "lon_a", nlon_atm, lon_dim_atm_id)
      call nc_error_handler(fid, ncstat, 10, "Defining lon dim atm")

      ncstat = nf90_def_dim(fid, "lat_a", nlat_atm, lat_dim_atm_id)
      call nc_error_handler(fid, ncstat, 11, "Defining lat dim atm")

      ncstat = nf90_def_dim(fid, "nlon_o", nlon_ocn, lon_dim_ocn_id)
      call nc_error_handler(fid, ncstat, 12, "Defining lon dim ocn")

      ncstat = nf90_def_dim(fid, "nlat_o", nlat_ocn, lat_dim_ocn_id)
      call nc_error_handler(fid, ncstat, 13, "Defining lat dim ocn")

      ncstat = nf90_def_dim(fid, "time", NF90_UNLIMITED, time_dim_id)
      call nc_error_handler(fid, ncstat, 14, "Defining time dim")

      !--- Define time coordinate
      ncstat = nf90_def_var(fid, "time", NF90_DOUBLE, (/time_dim_id/), time_id)
      call nc_error_handler(fid, ncstat, 60, "Defining time")
      ncstat = nf90_put_att( fid, time_id, "long_name", "Time" )
      call nc_error_handler(fid, ncstat, 60, "Defining time long_name")

      !--- Use reference time rather than start time to determine time origin
      !--- Start time will be 1 coupling interval before the start of the month but
      !--- the reference time will be the start of the coupling cycle (ie start of month)
      time_string = " "
      call ESMF_TimeGet(cpl_ref_time, timeStringISOFrac=time_string)
      do idx1=1,len_trim(time_string)
        if ( time_string(idx1:idx1) .eq. "T" ) time_string(idx1:idx1) = " "
      enddo
      ncstat = nf90_put_att( fid, time_id, "units", "days since "//trim(time_string) )
      call nc_error_handler(fid, ncstat, 60, "Defining time units")

      !--- Define ATM coordinate variables
      lon_atm_id = -1
      lat_atm_id = -1
      if ( len_trim(atm_grid_desc) > 0 ) then
        !--- Because the lon dimension size, nlon_atm, will always contain the overlap longitude
        !--- we must supply overlap = T in the call to read_scrip_grid_file so that the fields
        !--- mid_lon_atm, mid_lat_atm and mask_atm will be dimensioned consistent with nlon_atm
        grid_file_type = grid_desc_file_type(atm_grid_desc)
        select case (grid_file_type)
          case (1) !--- A SCRIP format file
            call read_scrip_grid_file(atm_grid_desc, grid_size=grid_size_atm, &
                   grid_dims=grid_dims_atm, mid_lon=mid_lon_atm, mid_lat=mid_lat_atm, &
                   imask=mask_atm, overlap=.true.)
          case (2) !--- A NEMO "mesh_mask" format file
            call read_NEMO_mesh_mask_file(atm_grid_desc, grid_size=grid_size_atm, &
                   grid_dims=grid_dims_atm, mid_lon=mid_lon_atm, mid_lat=mid_lat_atm, &
                   imask=mask_atm, overlap=.true.)
          case default
            write(6,*)'history_open: ',trim(atm_grid_desc), &
                '  Invalid grid_file_type = ',grid_file_type
            call xit("history_open",-1)
        end select
        if ( verbose > 0 ) then
          write(6,*)"history_open: ATM grid_size, grid_dims = ",grid_size_atm,grid_dims_atm
        endif
        if ( associated(mid_lon_atm) ) then
          !--- Define the ATM lon coordinate as a 1d array
          ncstat = nf90_def_var(fid, "lon_a", NF90_DOUBLE, &
                                (/lon_dim_atm_id/), lon_atm_id)
          call nc_error_handler(fid, ncstat, 8, "defining lon_a")
          ncstat = nf90_put_att( fid, lon_atm_id, "long_name", "Longitude (atm grid)" )
          call nc_error_handler(fid, ncstat, 8, "Defining lon_a long_name")
          ncstat = nf90_put_att( fid, lon_atm_id, "units", "degrees_east" )
          call nc_error_handler(fid, ncstat, 8, "Defining lon_a units")
          if (verbose > 1) then
            write(6,*)"history_open: nlon_atm, nlat_atm, shape(mid_lon_atm) = ", &
                      nlon_atm,nlat_atm,shape(mid_lon_atm)
          endif
        endif
        if ( associated(mid_lat_atm) ) then
          !--- Define the ATM lat coordinate as a 1d array
          ncstat = nf90_def_var(fid, "lat_a", NF90_DOUBLE, &
                                (/lat_dim_atm_id/), lat_atm_id)
          call nc_error_handler(fid, ncstat, 9, "defining lat_a")
          ncstat = nf90_put_att( fid, lat_atm_id, "long_name", "Latitude (atm grid)" )
          call nc_error_handler(fid, ncstat, 8, "Defining lat_a long_name")
          ncstat = nf90_put_att( fid, lat_atm_id, "units", "degrees_north" )
          call nc_error_handler(fid, ncstat, 8, "Defining lat_a units")
          if (verbose > 1) then
            write(6,*)"history_open: nlon_atm, nlat_atm, shape(mid_lat_atm) = ", &
                      nlon_atm,nlat_atm,shape(mid_lat_atm)
          endif
        endif

      else
        !--- Extract atm lon from the internal atm grid definition
        call get_atm_lon(mid_lon_atm, overlap=.true.)
        if ( size(mid_lon_atm,dim=1) /= nlon_atm .or. &
             size(mid_lon_atm,dim=2) /= nlat_atm ) then
          write(6,*)"history_open: Size mismatch for atm lon"
          write(6,*)"  expected ",nlon_atm,nlat_atm
          write(6,*)"     found ",shape(mid_lon_atm)
          call xit("history_open",-2)
        endif

        !--- Define the ATM lon coordinate as a 1d array
        ncstat = nf90_def_var(fid, "lon_a", NF90_DOUBLE, &
                              (/lon_dim_atm_id/), lon_atm_id)
        call nc_error_handler(fid, ncstat, 8, "defining lon_a")
        ncstat = nf90_put_att( fid, lon_atm_id, "long_name", "Longitude (atm grid)" )
        call nc_error_handler(fid, ncstat, 8, "Defining lon_a long_name")
        ncstat = nf90_put_att( fid, lon_atm_id, "units", "degrees_east" )
        call nc_error_handler(fid, ncstat, 8, "Defining lon_a units")
        if (verbose > -1) then
          write(6,*)"history_open: nlon_atm, nlat_atm, shape(mid_lon_atm) = ", &
                    nlon_atm,nlat_atm,shape(mid_lon_atm)
        endif

        !--- Extract atm lat from the internal atm grid definition
        call get_atm_lat(mid_lat_atm, overlap=.true.)
        if ( size(mid_lat_atm,dim=1) /= nlon_atm .or. &
             size(mid_lat_atm,dim=2) /= nlat_atm ) then
          write(6,*)"history_open: Size mismatch for atm lat"
          write(6,*)"  expected ",nlon_atm,nlat_atm
          write(6,*)"     found ",shape(mid_lat_atm)
          call xit("history_open",-2)
        endif

        !--- Define the ATM lat coordinate as a 1d array
        ncstat = nf90_def_var(fid, "lat_a", NF90_DOUBLE, &
                              (/lat_dim_atm_id/), lat_atm_id)
        call nc_error_handler(fid, ncstat, 9, "defining lat_a")
        ncstat = nf90_put_att( fid, lat_atm_id, "long_name", "Latitude (atm grid)" )
        call nc_error_handler(fid, ncstat, 8, "Defining lat_a long_name")
        ncstat = nf90_put_att( fid, lat_atm_id, "units", "degrees_north" )
        call nc_error_handler(fid, ncstat, 8, "Defining lat_a units")
        if (verbose > 1) then
          write(6,*)"history_open: nlon_atm, nlat_atm, shape(mid_lat_atm) = ", &
                    nlon_atm,nlat_atm,shape(mid_lat_atm)
        endif
      endif

      !--- Define OCN coordinate variables
      lon_ocn_id = -1
      lat_ocn_id = -1
      mask_ocn_id = -1
      if ( len_trim(ocn_grid_desc) > 0 ) then
        !--- Because the lon dimension size, nlon_ocn, will always contain the overlap longitude
        !--- we must supply overlap = T in the call to read_scrip_grid_file so that the fields
        !--- mid_lon_ocn, mid_lat_ocn and mask_ocn will be dimensioned consistent with nlon_ocn
        grid_file_type = grid_desc_file_type(ocn_grid_desc)
        select case (grid_file_type)
          case (1) !--- A SCRIP format file
            call read_scrip_grid_file(ocn_grid_desc, grid_size=grid_size_ocn, &
                   grid_dims=grid_dims_ocn, mid_lon=mid_lon_ocn, mid_lat=mid_lat_ocn, &
                   imask=mask_ocn, overlap=.true.)
          case (2) !--- A NEMO "mesh_mask" format file
            call read_NEMO_mesh_mask_file(ocn_grid_desc, grid_size=grid_size_ocn, &
                   grid_dims=grid_dims_ocn, mid_lon=mid_lon_ocn, mid_lat=mid_lat_ocn, &
                   imask=mask_ocn, overlap=.true.)
          case default
            write(6,*)'history_open: ',trim(ocn_grid_desc), &
                '  Invalid grid_file_type = ',grid_file_type
            call xit("history_open",-2)
        end select
        if ( verbose > 0 ) then
          write(6,*)"history_open: OCN grid_size, grid_dims = ",grid_size_ocn,grid_dims_ocn
        endif
        if ( associated(mid_lon_ocn) ) then
          !--- Define the OCN lon coordinate as a 2d array
          ncstat = nf90_def_var(fid, "lon_ocn", NF90_DOUBLE, &
                                (/lon_dim_ocn_id, lat_dim_ocn_id/), lon_ocn_id)
          call nc_error_handler(fid, ncstat, 10, "defining lon_ocn")
          ncstat = nf90_put_att( fid, lon_ocn_id, "long_name", "Longitude (ocn grid)" )
          call nc_error_handler(fid, ncstat, 10, "Defining lon_ocn long_name")
          ncstat = nf90_put_att( fid, lon_ocn_id, "units", "degrees_east" )
          call nc_error_handler(fid, ncstat, 10, "Defining lon_ocn units")
          if (verbose > 1) then
            write(6,*)"history_open: nlon_ocn, nlat_ocn, shape(mid_lon_ocn) = ", &
                      nlon_ocn,nlat_ocn,shape(mid_lon_ocn)
          endif
        endif
        if ( associated(mid_lat_ocn) ) then
          !--- Define the OCN lat coordinate as a 2d array
          ncstat = nf90_def_var(fid, "lat_ocn", NF90_DOUBLE, &
                                (/lon_dim_ocn_id, lat_dim_ocn_id/), lat_ocn_id)
          call nc_error_handler(fid, ncstat, 11, "defining lat_ocn")
          ncstat = nf90_put_att( fid, lat_ocn_id, "long_name", "Latitude (ocn grid)" )
          call nc_error_handler(fid, ncstat, 11, "Defining lat_ocn long_name")
          ncstat = nf90_put_att( fid, lat_ocn_id, "units", "degrees_north" )
          call nc_error_handler(fid, ncstat, 11, "Defining lat_ocn units")
          if (verbose > 1) then
            write(6,*)"history_open: nlon_ocn, nlat_ocn, shape(mid_lat_ocn) = ", &
                      nlon_ocn,nlat_ocn,shape(mid_lat_ocn)
          endif
        endif
        if ( associated(mask_ocn) ) then
          !--- Define the OCN mask as a 2d integer array
          ncstat = nf90_def_var(fid, "mask_ocn", NF90_INT, &
                                (/lon_dim_ocn_id, lat_dim_ocn_id/), mask_ocn_id)
          call nc_error_handler(fid, ncstat, 11, "defining mask_ocn")
          ncstat = nf90_put_att( fid, mask_ocn_id, "long_name", "Ocean Mask (1=ocean, 0=land)" )
          call nc_error_handler(fid, ncstat, 60, "Defining mask_ocn long_name")
          if (verbose > 1) then
            write(6,*)"history_open: nlon_ocn, nlat_ocn, shape(mask_ocn) = ", &
                      nlon_ocn,nlat_ocn,shape(mask_ocn)
          endif
        endif

      else
        !--- Extract ocn lon from the internal ocn grid definition
        call get_ocn_lon(mid_lon_ocn, overlap=.true.)
        if ( size(mid_lon_ocn,dim=1) /= nlon_ocn .or. &
             size(mid_lon_ocn,dim=2) /= nlat_ocn ) then
          write(6,*)"history_open: Size mismatch for ocn lon"
          write(6,*)"  expected ",nlon_ocn,nlat_ocn
          write(6,*)"     found ",shape(mid_lon_ocn)
          call xit("history_open",-2)
        endif

        !--- Define the OCN lon coordinate as a 2d array
        ncstat = nf90_def_var(fid, "lon_ocn", NF90_DOUBLE, &
                              (/lon_dim_ocn_id, lat_dim_ocn_id/), lon_ocn_id)
        call nc_error_handler(fid, ncstat, 10, "defining lon_ocn")
        ncstat = nf90_put_att( fid, lon_ocn_id, "long_name", "Longitude (ocn grid)" )
        call nc_error_handler(fid, ncstat, 10, "Defining lon_ocn long_name")
        ncstat = nf90_put_att( fid, lon_ocn_id, "units", "degrees_east" )
        call nc_error_handler(fid, ncstat, 10, "Defining lon_ocn units")
        if (verbose > 1) then
          write(6,*)"history_open: nlon_ocn, nlat_ocn, shape(mid_lon_ocn) = ", &
                    nlon_ocn,nlat_ocn,shape(mid_lon_ocn)
        endif

        !--- Extract ocn lat from the internal ocn grid definition
        call get_ocn_lat(mid_lat_ocn, overlap=.true.)
        if ( size(mid_lat_ocn,dim=1) /= nlon_ocn .or. &
             size(mid_lat_ocn,dim=2) /= nlat_ocn ) then
          write(6,*)"history_open: Size mismatch for ocn lat"
          write(6,*)"  expected ",nlon_ocn,nlat_ocn
          write(6,*)"     found ",shape(mid_lat_ocn)
          call xit("history_open",-2)
        endif

        !--- Define the OCN lat coordinate as a 2d array
        ncstat = nf90_def_var(fid, "lat_ocn", NF90_DOUBLE, &
                              (/lon_dim_ocn_id, lat_dim_ocn_id/), lat_ocn_id)
        call nc_error_handler(fid, ncstat, 11, "defining lat_ocn")
        ncstat = nf90_put_att( fid, lat_ocn_id, "long_name", "Latitude (ocn grid)" )
        call nc_error_handler(fid, ncstat, 11, "Defining lat_ocn long_name")
        ncstat = nf90_put_att( fid, lat_ocn_id, "units", "degrees_north" )
        call nc_error_handler(fid, ncstat, 11, "Defining lat_ocn units")
        if (verbose > 1) then
          write(6,*)"history_open: nlon_ocn, nlat_ocn, shape(mid_lat_ocn) = ", &
                    nlon_ocn,nlat_ocn,shape(mid_lat_ocn)
        endif

        !--- Extract ocn mask from the internal ocn grid definition
        call get_ocn_mask(mask_ocn, overlap=.true.)
        if ( size(mask_ocn,dim=1) /= nlon_ocn .or. &
             size(mask_ocn,dim=2) /= nlat_ocn ) then
          write(6,*)"history_open: Size mismatch for ocn mask"
          write(6,*)"  expected ",nlon_ocn,nlat_ocn
          write(6,*)"     found ",shape(mask_ocn)
          call xit("history_open",-2)
        endif

        !--- Define the OCN mask as a 2d integer array
        ncstat = nf90_def_var(fid, "mask_ocn", NF90_INT, &
                              (/lon_dim_ocn_id, lat_dim_ocn_id/), mask_ocn_id)
        call nc_error_handler(fid, ncstat, 11, "defining mask_ocn")
        ncstat = nf90_put_att( fid, mask_ocn_id, "long_name", "Ocean Mask (1=ocean, 0=land)" )
        call nc_error_handler(fid, ncstat, 60, "Defining mask_ocn long_name")
        if (verbose > 1) then
          write(6,*)"history_open: nlon_ocn, nlat_ocn, shape(mask_ocn) = ", &
                    nlon_ocn,nlat_ocn,shape(mask_ocn)
        endif

      endif

      !--- Define the AGCM fractional land mask
      atm_fland_id = -1
      if ( associated(atm_fland) ) then
        !--- Define the ATM mask as a 2d real*8 array
        if ( size(atm_fland) /= nlon_atm*nlat_atm ) then
          write(6,*)"history_open: shape missmatch for atm_fland"
          write(6,*)"    atm_fland nx*ny: ",size(atm_fland)
          write(6,*)"    expected  nx*ny: ",nlon_atm*nlat_atm
          call xit("history_open",-4)
        endif
        ncstat = nf90_def_var(fid, "fland_atm", NF90_DOUBLE, &
                              (/lon_dim_atm_id, lat_dim_atm_id/), atm_fland_id)
        call nc_error_handler(fid, ncstat, 9, "defining atm_fland")
        ncstat = nf90_put_att( fid, atm_fland_id, "long_name", "AGCM Fractional Land Mask" )
        call nc_error_handler(fid, ncstat, 60, "Defining fland_atm long_name")
        ncstat = nf90_put_att( fid, atm_fland_id, "units", "1" )
        call nc_error_handler(fid, ncstat, 61, "Defining fland_atm units")
        if (verbose > 1) then
          write(6,*)"history_open: nlon_atm, nlat_atm, shape(atm_fland) = ", &
                    nlon_atm,nlat_atm,shape(atm_fland)
        endif
      endif

      !--- Define the AGCM grid cell area
      atm_grid_cell_area_id = -1
      if ( associated(atm_grid_cell_area) ) then
        !--- Define the ATM mask as a 2d real*8 array
        !--- Allow atm_grid_cell_area to contain or not contain the overlap lon
        !--- The field written to the history file will always have shape (nlon_atm,nlat_atm)
        found = (size(atm_grid_cell_area,dim=1) == nlon_atm .or. &
                 size(atm_grid_cell_area,dim=1) == nlon_atm-1) .and. &
                 size(atm_grid_cell_area,dim=2) == nlat_atm
        if ( .not. found ) then
          write(6,*)"history_open: shape missmatch for atm_grid_cell_area"
          write(6,*)"    atm_grid_cell_area nx,ny: ",size(atm_grid_cell_area,dim=1),size(atm_grid_cell_area,dim=2)
          write(6,*)"    expected           nx,ny: ",nlon_atm,nlat_atm," or ",nlon_atm-1,nlat_atm
          call xit("history_open",-4)
        endif
        ncstat = nf90_def_var(fid, "atm_grid_cell_area", NF90_DOUBLE, &
                              (/lon_dim_atm_id, lat_dim_atm_id/), atm_grid_cell_area_id)
        call nc_error_handler(fid, ncstat, 9, "defining atm_grid_cell_area")
        ncstat = nf90_put_att( fid, atm_grid_cell_area_id, "long_name", "AGCM Grid Cell Area" )
        call nc_error_handler(fid, ncstat, 60, "Defining atm_grid_cell_area long_name")
        if ( maxval(atm_grid_cell_area) > 1e4 ) then
          !--- Assume units are meters squared
          ncstat = nf90_put_att( fid, atm_grid_cell_area_id, "units", "m2" )
        else
          !--- Assume units are dimensionless
          ncstat = nf90_put_att( fid, atm_grid_cell_area_id, "units", "1" )
        endif
        call nc_error_handler(fid, ncstat, 61, "Defining atm_grid_cell_area units")
        if (verbose > 1) then
          write(6,*)"history_open: nlon_atm, nlat_atm, shape(atm_grid_cell_area) = ", &
                    nlon_atm,nlat_atm,shape(atm_grid_cell_area)
        endif
      endif

      !--- Define the ocean grid cell area
      ocn_grid_cell_area_id = -1
      if ( associated(ocn_grid_cell_area) ) then
        !--- Define the OCN mask as a 2d real*8 array
        !--- Allow ocn_grid_cell_area to contain or not contain the overlap lons
        !--- The field written to the history file will always have shape (nlon_ocn,nlat_ocn)
        found = (size(ocn_grid_cell_area,dim=1) == nlon_ocn .or. &
                 size(ocn_grid_cell_area,dim=1) == nlon_ocn-2) .and. &
                 size(ocn_grid_cell_area,dim=2) == nlat_ocn
        if ( .not. found ) then
          write(6,*)"history_open: shape missmatch for ocn_grid_cell_area"
          write(6,*)"    ocn_grid_cell_area nx,ny: ",size(ocn_grid_cell_area,dim=1),size(ocn_grid_cell_area,dim=2)
          write(6,*)"    expected           nx,ny: ",nlon_ocn,nlat_ocn," or ",nlon_ocn-2,nlat_ocn
          call xit("history_open",-4)
        endif
        ncstat = nf90_def_var(fid, "ocn_grid_cell_area", NF90_DOUBLE, &
                              (/lon_dim_ocn_id, lat_dim_ocn_id/), ocn_grid_cell_area_id)
        call nc_error_handler(fid, ncstat, 9, "defining ocn_grid_cell_area")
        ncstat = nf90_put_att( fid, ocn_grid_cell_area_id, "long_name", "NEMO Grid Cell Area" )
        call nc_error_handler(fid, ncstat, 60, "Defining ocn_grid_cell_area long_name")
        if ( maxval(ocn_grid_cell_area) > 1e4 ) then
          !--- Assume units are meters squared
          ncstat = nf90_put_att( fid, ocn_grid_cell_area_id, "units", "m2" )
        else
          !--- Assume units are dimensionless
          ncstat = nf90_put_att( fid, ocn_grid_cell_area_id, "units", "1" )
        endif
        call nc_error_handler(fid, ncstat, 61, "Defining ocn_grid_cell_area units")
        if (verbose > 1) then
          write(6,*)"history_open: nlon_ocn, nlat_ocn, shape(ocn_grid_cell_area) = ", &
                    nlon_ocn,nlat_ocn,shape(ocn_grid_cell_area)
        endif
      endif

      !--- All other history variables will be defined the first time
      !--- they are written to the history file

      !--- Get out of define mode
      ncstat = nf90_enddef(fid)
      call nc_error_handler(fid, ncstat, 14)

      !--- Write coordinate values to the history file
      if ( lon_atm_id > -1 ) then
        if ( .not. associated(mid_lon_atm) ) then
          write(6,'(3a)')"history_open: Data values are missing for lon_atm."
          call xit("history_open",-3)
        endif
        ncstat = nf90_put_var( fid, lon_atm_id, mid_lon_atm(1:nlon_atm,1) )
        call nc_error_handler(fid, ncstat, 51, "assigning lon_atm")
      endif

      if ( lat_atm_id > -1 ) then
        if ( .not. associated(mid_lat_atm) ) then
          write(6,'(3a)')"history_open: Data values are missing for lat_atm."
          call xit("history_open",-4)
        endif
        ncstat = nf90_put_var( fid, lat_atm_id, mid_lat_atm(1,1:nlat_atm) )
        call nc_error_handler(fid, ncstat, 52, "assigning lat_atm")
      endif

      if ( lon_ocn_id > -1 ) then
        if ( .not. associated(mid_lon_ocn) ) then
          write(6,'(3a)')"history_open: Data values are missing for lon_ocn."
          call xit("history_open",-5)
        endif
        ncstat = nf90_put_var( fid, lon_ocn_id, mid_lon_ocn )
        call nc_error_handler(fid, ncstat, 53, "assigning lon_ocn")
      endif

      if ( lat_ocn_id > -1 ) then
        if ( .not. associated(mid_lat_ocn) ) then
          write(6,'(3a)')"history_open: Data values are missing for lat_ocn."
          call xit("history_open",-6)
        endif
        ncstat = nf90_put_var( fid, lat_ocn_id, mid_lat_ocn )
        call nc_error_handler(fid, ncstat, 54, "assigning lat_ocn")
      endif

      if ( mask_ocn_id > -1 ) then
        if ( .not. associated(mask_ocn) ) then
          write(6,'(3a)')"history_open: Data values are missing for mask_ocn."
          call xit("history_open",-8)
        endif
        ncstat = nf90_put_var( fid, mask_ocn_id, mask_ocn )
        call nc_error_handler(fid, ncstat, 51, "assigning mask_ocn")

        !--- Define a cpl_vinfo_t variable to use for the OCN mask
        allocate( tmp_vinfo )
        tmp_vinfo%name     = "Mask_ocn"
        tmp_vinfo%grid     = "ocn"
        tmp_vinfo%tag      = cpl_var_new_tag(5000)
        tmp_vinfo%size     = size(mask_ocn, dim=1) * size(mask_ocn, dim=2)
        tmp_vinfo%nlon     = size(mask_ocn, dim=1) - olap_o
        tmp_vinfo%nlat     = size(mask_ocn, dim=2)
        tmp_vinfo%olap     = olap_o
        tmp_vinfo%ncopy    = 1
        tmp_vinfo%freq     = 0
        tmp_vinfo%ccc_name = "OMSK"
        tmp_vinfo%regrid%method   = 'bilinear'
        tmp_vinfo%regrid%norm_opt = 'fracarea'
        tmp_vinfo%regrid%masked   = .false.
        !--- Add this variable defintion
        call add_cpl_var(tmp_vinfo, exists=0)
        nullify( tmp_vinfo )

        !--- Add the OCN mask to internal storage
        call store_AppendRecord("Mask_ocn", ocn_gid, mask_ocn, create=.true.)
      endif

      !--- Write agcm fractional land mask to the history file
      if ( atm_fland_id > -1 ) then
        if ( .not. associated(atm_fland) ) then
          write(6,'(3a)')"history_open: Data values are missing for atm_fland."
          call xit("history_open",-7)
        endif
        if ( associated(wrk2d) ) deallocate(wrk2d)
        allocate( wrk2d(nlon_atm,nlat_atm) )
        !--- Assign to a 2d array for the write to the history
        nwrds = nlon_atm*nlat_atm
        wrk2d = reshape( atm_fland(1:nwrds), (/nlon_atm, nlat_atm/) )
        ncstat = nf90_put_var( fid, atm_fland_id, wrk2d )
        call nc_error_handler(fid, ncstat, 51, "assigning atm_fland")
        deallocate(wrk2d)
      endif

      !--- Write agcm grid cell area to the history file
      if ( atm_grid_cell_area_id > -1 ) then
        if ( associated(atm_grid_cell_area) ) then
          if ( size(atm_grid_cell_area,dim=1) == nlon_atm .and. &
               size(atm_grid_cell_area,dim=2) == nlat_atm) then
            ncstat = nf90_put_var( fid, atm_grid_cell_area_id, atm_grid_cell_area )
            call nc_error_handler(fid, ncstat, 51, "assigning atm_grid_cell_area")
          else if ( size(atm_grid_cell_area,dim=1) == nlon_atm-1 .and. &
                    size(atm_grid_cell_area,dim=2) == nlat_atm) then
            !--- Add the overlap longitude to the output array
            if ( associated(wrk2d) ) deallocate(wrk2d)
            allocate( wrk2d(nlon_atm,nlat_atm) )
            wrk2d(1:nlon_atm-1,:) = atm_grid_cell_area(1:nlon_atm-1,:)
            wrk2d(nlon_atm,:) = atm_grid_cell_area(1,:)
            ncstat = nf90_put_var( fid, atm_grid_cell_area_id, wrk2d )
            call nc_error_handler(fid, ncstat, 51, "assigning atm_grid_cell_area")
            deallocate(wrk2d)
          endif
        else
          write(6,'(3a)')"history_open: Data values are missing for atm_grid_cell_area."
          call xit("history_open",-7)
        endif
      endif

      !--- Write ocean grid cell area to the history file
      if ( ocn_grid_cell_area_id > -1 ) then
        if ( associated(ocn_grid_cell_area) ) then
          if ( size(ocn_grid_cell_area,dim=1) == nlon_ocn .and. &
               size(ocn_grid_cell_area,dim=2) == nlat_ocn) then
            ncstat = nf90_put_var( fid, ocn_grid_cell_area_id, ocn_grid_cell_area )
            call nc_error_handler(fid, ncstat, 51, "assigning ocn_grid_cell_area")
          else if ( size(ocn_grid_cell_area,dim=1) == nlon_ocn-2 .and. &
                    size(ocn_grid_cell_area,dim=2) == nlat_ocn) then
            !--- Add the overlap longitude to the output array
            if ( associated(wrk2d) ) deallocate(wrk2d)
            allocate( wrk2d(nlon_ocn,nlat_ocn) )
            wrk2d(1:nlon_ocn-2,:) = ocn_grid_cell_area(1:nlon_ocn-2,:)
            wrk2d(nlon_ocn-1,:)   = ocn_grid_cell_area(1,:)
            wrk2d(nlon_ocn,:)     = ocn_grid_cell_area(2,:)
            ncstat = nf90_put_var( fid, ocn_grid_cell_area_id, wrk2d )
            call nc_error_handler(fid, ncstat, 51, "assigning ocn_grid_cell_area")
            deallocate(wrk2d)
          endif
        else
          write(6,'(3a)')"history_open: Data values are missing for ocn_grid_cell_area."
          call xit("history_open",-7)
        endif
      endif

      if ( verbose > -1 ) then
        write(6,'(2a)')"history_open: Create coupler history file ",trim(hist_file)
        call flush(6)
      endif

    end subroutine history_open

end module history

