module restart
  !--- Restart related procedures for CanCPL

  use com_cpl
  use ccc_esmf
  use ccc_netcdf
  use store
  use ccc_io

  implicit none

  private
  public :: restart_read
  public :: restart_write
  public :: cpl_rs_var_name_atm
  public :: cpl_rs_var_name_ocn

  character(len=256), public :: rs_file_in  = " "
  character(len=256), public :: rs_file_out = " "

  !--- Define grid sizes used in the restart file
  !--- 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

  !--- Define the list of variable to be saved in the restart file
  !--- Pad all initialized character variables with blanks to avoid a quirk of certain
  !--- compilers (e.g. xlf) wherein all initialized elements of a character array are
  !--- truncated to the length of the first element in the array

  !--- Variables on the ATM grid
  integer, save :: nvar_atm = 0
  character(32), pointer, save :: cpl_rs_var_name_atm(:) => null()
  integer, pointer, save       :: cpl_rs_var_gid_atm(:)  => null()

  !--- Variables on the OCN grid
  integer, save :: nvar_ocn = 0
  character(32), pointer, save :: cpl_rs_var_name_ocn(:) => null()
  integer, pointer, save       :: cpl_rs_var_gid_ocn(:)  => null()

  character(len=256) :: strng = " "

  contains

    !***************************************************************************
    !--- Initialize cpl_rs_var_name_atm and cpl_rs_var_name_ocn
    !***************************************************************************
    subroutine restart_init()
      integer :: idx1, idx2, idx3, verbose=1
      character(32) :: tmp_names(200)
      integer :: tmp_gids(200), gid
      logical :: found

      if (nvar_ocn == 0) then
        !--- Define the list of variables on the ocn grid

        if ( .not. associated(event) ) then
          write(6,*)"restart_init: event is not defined."
          write(6,*)"    You must define cpl_rs_var_name_atm and cpl_rs_var_name_ocn manually."
          call xit("restart_init",-1)
        endif

        !--- Initialize a temporary list
        tmp_names(:) = " "
        tmp_gids(:) = 0

        !--- Extract ocn variable names from the event list
        do idx1=1,nevents
          if ( event(idx1)%is_msg ) then
            !--- Ignore message variables
            cycle
          endif

          !--- Count the number of variables to be sent to or received from the ocean
          !--- and populate tmp_names with their names
          if ( trim(adjustl(event(idx1)%put_model)) .eq. "ocn" .or. &
               trim(adjustl(event(idx1)%get_model)) .eq. "ocn" ) then
            if ( nvar_ocn == 0 ) then
              nvar_ocn = nvar_ocn + 1
              tmp_names(nvar_ocn) = trim(adjustl(event(idx1)%put_var))
              if ( trim(adjustl(event(idx1)%put_model)) .eq. "ocn" ) then
                !--- Sending data to ocean, assume data comes from the atm model
                tmp_gids(nvar_ocn) = atm_gid
              else
                !--- Sending data from the ocean
                tmp_gids(nvar_ocn) = ocn_gid
              endif
              !--- Reset the group ID for certain variables
              gid = tmp_gids(nvar_ocn)
              tmp_gids(nvar_ocn) = filter_group_ID(event(idx1)%put_var, gid)
            else
              !--- Make sure there are no duplicate names
              found = .false.
              do idx2=1,nvar_ocn
                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
                nvar_ocn = nvar_ocn + 1
                tmp_names(nvar_ocn) = trim(adjustl(event(idx1)%put_var))
                if ( trim(adjustl(event(idx1)%put_model)) .eq. "ocn" ) then
                  !--- Sending data to ocean, assume data comes from the atm model
                  tmp_gids(nvar_ocn) = atm_gid
                else
                  !--- Sending data from the ocean
                  tmp_gids(nvar_ocn) = ocn_gid
                endif
                !--- Reset the group ID for certain variables
                gid = tmp_gids(nvar_ocn)
                tmp_gids(nvar_ocn) = filter_group_ID(event(idx1)%put_var, gid)
              endif
            endif
          endif
        enddo

        if ( nvar_ocn > 0 ) then
          !--- Allocate and assign cpl_rs_var_name_ocn
          allocate( cpl_rs_var_name_ocn(nvar_ocn) )
          cpl_rs_var_name_ocn(1:nvar_ocn) = tmp_names(1:nvar_ocn)
          allocate( cpl_rs_var_gid_ocn(nvar_ocn) )
          cpl_rs_var_gid_ocn(1:nvar_ocn) = tmp_gids(1:nvar_ocn)
        else
          !--- Set nvar_ocn to -1 to indicate that there are no ocn variables
          !--- to be added to the restart file
          nvar_ocn = -1
        endif

        if ( verbose > 0 ) then
          write(6,'(a,i3)')'restart_init: nvar_ocn = ',nvar_ocn
          if ( nvar_ocn > 0 ) write(6,'(8(2x,a12," gid=",i2))') &
             (cpl_rs_var_name_ocn(idx1),cpl_rs_var_gid_ocn(idx1),idx1=1,nvar_ocn)
          call flush(6)
        endif

      endif

      if (nvar_atm == 0) then
        !--- Define the list of variables on the atm grid

        if ( .not. associated(event) ) then
          write(6,*)"restart_init: event is not defined."
          write(6,*)"    You must define cpl_rs_var_name_atm and cpl_rs_var_name_atm manually."
          call xit("restart_init",-1)
        endif

        !--- Initialize a temporary list
        tmp_names(:) = " "
        tmp_gids(:) = 0

        !--- Extract atm variable names from the event list
        do idx1=1,nevents
          if ( event(idx1)%is_msg ) then
            !--- Ignore message variables
            cycle
          endif

          !--- Count the number of variables to be sent to or received from the atm
          !--- and populate tmp_names with their names
          if ( trim(adjustl(event(idx1)%put_model)) .eq. "atm" .or. &
               trim(adjustl(event(idx1)%get_model)) .eq. "atm" ) then
            if ( nvar_atm == 0 ) then
              nvar_atm = nvar_atm + 1
              tmp_names(nvar_atm) = trim(adjustl(event(idx1)%put_var))
              !--- Determine the group ID associated with this variable
              if ( trim(adjustl(event(idx1)%put_model)) .eq. "atm" ) then
                !--- Sending data to agcm, assume data comes from the ocean model
                tmp_gids(nvar_atm) = ocn_gid
              else
                !--- Sending data from the agcm
                tmp_gids(nvar_atm) = atm_gid
              endif
              !--- Reset the group ID for certain variables
              gid = tmp_gids(nvar_atm)
              tmp_gids(nvar_atm) = filter_group_ID(event(idx1)%put_var, gid)
            else
              !--- Make sure there are no duplicate names
              found = .false.
              do idx2=1,nvar_atm
                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
                nvar_atm = nvar_atm + 1
                tmp_names(nvar_atm) = trim(adjustl(event(idx1)%put_var))
                !--- Determine the group ID associated with this variable
                if ( trim(adjustl(event(idx1)%put_model)) .eq. "atm" ) then
                  !--- Sending data to agcm, assume data comes from the ocean model
                  tmp_gids(nvar_atm) = ocn_gid
                else
                  !--- Sending data from the agcm
                  tmp_gids(nvar_atm) = atm_gid
                endif
                !--- Reset the group ID for certain variables
                gid = tmp_gids(nvar_atm)
                tmp_gids(nvar_atm) = filter_group_ID(event(idx1)%put_var, gid)
              endif
            endif
          endif
        enddo

        if ( nvar_atm > 0 ) then
          !--- Allocate and assign cpl_rs_var_name_atm
          allocate( cpl_rs_var_name_atm(nvar_atm) )
          cpl_rs_var_name_atm(1:nvar_atm) = tmp_names(1:nvar_atm)
          allocate( cpl_rs_var_gid_atm(nvar_atm) )
          cpl_rs_var_gid_atm(1:nvar_atm) = tmp_gids(1:nvar_atm)
        else
          !--- Set nvar_atm to -1 to indicate that there are no atm variables
          !--- to be added to the restart file
          nvar_atm = -1
        endif

        if ( verbose > 0 ) then
          write(6,'(a,i3)')'restart_init: nvar_atm = ',nvar_atm
          if ( nvar_atm > 0 ) write(6,'(8(2x,a12," gid=",i2))') &
             (cpl_rs_var_name_atm(idx1),cpl_rs_var_gid_atm(idx1),idx1=1,nvar_atm)
          call flush(6)
        endif

      endif

    end subroutine restart_init

    !***************************************************************************
    !--- Read a CanCPL restart file
    !***************************************************************************
    subroutine restart_read(rs_file, time_string, init)
      character(*), intent(in) :: rs_file
      character(*), intent(out) :: time_string
      integer, intent(in), optional :: init

      !--- Local
      integer(kind=4) :: rs_fid, ncstat, vid
      integer :: idx1, nwrds, gid
      integer :: local_init
      logical :: found
      character(32) :: var_name
      real(kind=8), pointer :: vals(:) => null()
      real(kind=8), pointer :: vals2d(:,:) => null()

      !--- Determine if restart values are initialized internally (and if so, how)
      !--- or read from the restart file
      if ( present(init) ) then
        !--- Set a flag to indicate the type of initialization
        local_init = init
      else
        !--- The default is to read from the restart file
        local_init = 0
      endif

      !--- Ensure the list of restart variables is defined
      call restart_init()

      !--- Open the restart file for reading
      ncstat = nf90_open(rs_file, NF90_NOWRITE, rs_fid)
      call nc_error_handler(rs_fid, ncstat, 1, "Attempting to read restart file "//trim(rs_file))

      !--- Read the global attribute restart_time into time_string
      strng = " "
      ncstat = nf90_get_att(rs_fid, NF90_GLOBAL, "restart_time", strng)
      call nc_error_handler(rs_fid, ncstat, 2, "reading global attribute restart_time")
      time_string = " "
      time_string = trim(strng)

      if ( nvar_atm > 0 ) then
        !--- Allocate space for 1D and 2D versions of a variable on the ATM grid
        allocate( vals(nlon_atm*nlat_atm) )
        allocate( vals2d(nlon_atm,nlat_atm) )

        !--- The total number of words in any ATM restart field
        nwrds = nlon_atm*nlat_atm

        !--- Read all restart variables that are on the ATM grid
        do idx1=1,size(cpl_rs_var_name_atm)
          var_name = " "
          var_name = trim(cpl_rs_var_name_atm(idx1))
          gid = cpl_rs_var_gid_atm(idx1)

          !--- Read the current variable from a NEMO restart file, if possible
          !--- For certain variables on the ATM grid, read_nemo_restart will read a
          !--- corresponding NEMO variable and interpolate it onto the ATM grid
          call read_nemo_restart(var_name, gid, found)

          !--- If the NEMO restart file was present and the variable was found in it
          !--- then ignore the remainder of this loop
          if ( found ) cycle

          !--- Read the current variable from the coupler restart file
          ncstat =  nf90_inq_varid(rs_fid, trim(cpl_rs_var_name_atm(idx1)), vid)
          if (ncstat /= NF90_NOERR) then
            !--- The current variable does not exists in the restart file
            if ( cpl_rs_abort_if_missing_field ) then
              write(6,*)"restart_read: The variable ",trim(cpl_rs_var_name_atm(idx1)), &
                " is missing from the coupler restart file."
              call xit("restart_read",-1)
            endif
            call restart_get_initial_values(var_name, nlon_atm, nlat_atm, vals2d)
            !--- Copy the values read from the netcdf file to a 1d array
            vals(1:nwrds) = reshape( vals2d(1:nlon_atm,1:nlat_atm), (/ nwrds /) )
            write(6,*)"restart_read: ",trim(cpl_rs_var_name_atm(idx1)), &
                " is not in the coupler restart file. Inserting default values."
          else
            !--- The variable exists so read it
            ncstat = nf90_get_var(rs_fid, vid, vals2d)
            call nc_error_handler(rs_fid, ncstat, 3)

            !--- Copy the values read from the netcdf file to a 1d array
            vals(1:nwrds) = reshape( vals2d(1:nlon_atm,1:nlat_atm), (/ nwrds /) )
          endif

          !!!!!!!!!!!!! The AGCM will die without this !!!!!!!!!!!!!!
          !--- The AGCM will die, typically via xit(-1) from PHYSICI
          !--- The zero values in GT cause a problem in SFCPROC when
          !--- it does its convoluted tiled surface temp calculation.
          !--- The problem arises in OIFPRP10 when it divides by the
          !--- internally calculated GTGAT (gathered from GTROT)
          !--- and the GTGAT array contains zeros, resulting in NaNs.
          if ( trim(adjustl(var_name)) .eq. "GT_atm" ) then
            where( vals < 10.0 ) vals = 273.15
          endif
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#undef DBG_bfb
#ifdef DBG_bfb
strng = " "
write(strng,'(2a)')"restart_var_",trim(cpl_rs_var_name_atm(idx1))
call dump_farray(trim(strng),vals(1:nwrds),idx=cpl_elapsed_time_secs)
#endif

          !--- Add this field to internal storage
          call store_AppendRecord(cpl_rs_var_name_atm(idx1), cpl_rs_var_gid_atm(idx1), &
                                  vals, create=.true.)

        enddo
        if ( associated(vals)   ) deallocate( vals )
        if ( associated(vals2d) ) deallocate( vals2d )
      endif

      if ( nvar_ocn > 0 ) then
        !--- Allocate space for 1D and 2D versions of a variable on the OCN grid
        allocate( vals(nlon_ocn*nlat_ocn) )
        allocate( vals2d(nlon_ocn,nlat_ocn) )

        !--- The total number of words in any OCN restart field
        nwrds = nlon_ocn*nlat_ocn

        !--- Read all restart variables that are on the OCN grid
        do idx1=1,size(cpl_rs_var_name_ocn)
          var_name = " "
          var_name = trim(cpl_rs_var_name_ocn(idx1))
          gid = cpl_rs_var_gid_ocn(idx1)

          !--- Read the current variable from a NEMO restart file, if possible
          call read_nemo_restart(var_name, gid, found)

          !--- If the NEMO restart file was present and the variable was found in it
          !--- then ignore the remainder of this loop
          if ( found ) cycle

          !--- Read the current variable from the restart file
          ncstat =  nf90_inq_varid(rs_fid, trim(cpl_rs_var_name_ocn(idx1)), vid)
          if (ncstat /= NF90_NOERR) then
            !--- The current variable does not exists in the restart file
            if ( cpl_rs_abort_if_missing_field ) then
              write(6,*)"restart_read: The variable ",trim(cpl_rs_var_name_ocn(idx1)), &
                " is missing from the coupler restart file."
              call xit("restart_read",-1)
            endif
            call restart_get_initial_values(var_name, nlon_ocn, nlat_ocn, vals2d)
            !--- Copy the values read from the netcdf file to a 1d array
            vals(1:nwrds) = reshape( vals2d(1:nlon_ocn,1:nlat_ocn), (/ nwrds /) )
            write(6,*)"restart_read: ",trim(cpl_rs_var_name_ocn(idx1)), &
                " is not in the coupler restart file. Inserting default values."
          else
            !--- The variable exists so read it
            ncstat = nf90_get_var(rs_fid, vid, vals2d)
            call nc_error_handler(rs_fid, ncstat, 3)

            !--- Copy the values read from the netcdf file to a 1d array
            vals(1:nwrds) = reshape( vals2d(1:nlon_ocn,1:nlat_ocn), (/ nwrds /) )
          endif

#ifdef DBG_bfb
strng = " "
write(strng,'(2a)')"restart_var_",trim(cpl_rs_var_name_ocn(idx1))
call dump_farray(trim(strng),vals(1:nwrds),idx=cpl_elapsed_time_secs)
#endif

          !--- Add this field to internal storage
          call store_AppendRecord(cpl_rs_var_name_ocn(idx1), cpl_rs_var_gid_ocn(idx1), &
                                  vals, create=.true.)

        enddo
        if ( associated(vals)   ) deallocate( vals )
        if ( associated(vals2d) ) deallocate( vals2d )
      endif

      !--- Close the netcdf file
      ncstat = nf90_close(rs_fid)
      call nc_error_handler(rs_fid, ncstat, 99)

    end subroutine restart_read

    !***************************************************************************
    !--- Write a CanCPL restart file
    !***************************************************************************
    subroutine restart_write(rs_file, time_string, initial, save_to_datapath)
      character(*), intent(in) :: rs_file
      character(*), intent(in) :: time_string
      logical, intent(in), optional :: initial
      logical, intent(in), optional :: save_to_datapath

      !--- Local
      integer(kind=4) :: rs_fid, ncstat, vid
      integer :: idx1, nwrds
      character(32) :: var_name
      character(128) :: cs_file
      logical :: write_initial_values
      !--- Space for data that is read from storage and written to the restart
      real(kind=8), pointer :: vals(:) => null()
      real(kind=8), pointer :: vals2d(:,:) => null()
      integer(kind=8) :: time_out

      !--- Determine if initial values are written to the restart file
      !--- rather than field values read from internal storage
      if ( present(initial) ) then
        write_initial_values = initial
      else
        write_initial_values = .false.
      endif

      !--- Create a restart file, define all variables and return the netcdf file id
      call restart_create(rs_file, rs_fid, time_string)

      if ( nvar_atm > 0 ) then
        !--- Write atm fields to the restart file
        allocate( vals2d(nlon_atm,nlat_atm) )
        do idx1=1,size(cpl_rs_var_name_atm)
          var_name = trim(adjustl(cpl_rs_var_name_atm(idx1)))
          time_out = 0

          if ( write_initial_values ) then
            !--- Use initial values for restart fields
            call restart_get_initial_values(var_name, nlon_atm, nlat_atm, vals2d)

          else
            !--- Get field data from internal storage
            !--- Get the last record that was written to storage for the current name
            nwrds = store_readRecord(var_name, cpl_rs_var_gid_atm(idx1), vals, idx=-1, time_out=time_out)

            if ( nwrds /= nlon_atm*nlat_atm ) then
              write(6,'(3a)')"restart_write: Size mismatch reading ", &
                             trim(cpl_rs_var_name_atm(idx1))," from internal storage."
              write(6,*)"   Expected ",nlon_atm*nlat_atm,"  Found ",nwrds
              call xit("restart_write",-1)
            endif

            !--- Assign these values to a 2d array for the write to the restart
            vals2d(:,:) = reshape( vals(1:nwrds), (/nlon_atm,nlat_atm/) )
          endif

          !--- Write this field to the restart file
          ncstat =  nf90_inq_varid(rs_fid, trim(cpl_rs_var_name_atm(idx1)), vid)
          call nc_error_handler(rs_fid, ncstat, 1)
          ncstat = nf90_put_var( rs_fid, vid, vals2d )
          call nc_error_handler(rs_fid, ncstat, 2, "writing "//trim(var_name))

          write(6,'(3a,i10,2a)')"restart_write: Wrote ",trim(cpl_rs_var_name_atm(idx1)), &
                         " at time ",time_out," to ",trim(rs_file)
          call flush(6)
        enddo
        if ( associated(vals)   ) deallocate( vals )
        if ( associated(vals2d) ) deallocate( vals2d )
      endif

      if ( nvar_ocn > 0 ) then
        !--- Write ocn fields to the restart file
        allocate( vals2d(nlon_ocn,nlat_ocn) )
        do idx1=1,size(cpl_rs_var_name_ocn)
          var_name = trim(adjustl(cpl_rs_var_name_ocn(idx1)))
          time_out = 0

          if ( write_initial_values ) then
            !--- Use initial values for restart fields
            call restart_get_initial_values(var_name, nlon_ocn, nlat_ocn, vals2d)

          else
            !--- Get field data from internal storage
            !--- Get the last record that was written to storage for the current name
            nwrds = store_readRecord(var_name, cpl_rs_var_gid_ocn(idx1), vals, idx=-1, time_out=time_out)

            if ( nwrds /= nlon_ocn*nlat_ocn ) then
              write(6,'(3a)')"restart_write: Size mismatch reading ", &
                             trim(cpl_rs_var_name_ocn(idx1))," from internal storage."
              write(6,*)"   Expected ",nlon_ocn*nlat_ocn,"  Found ",nwrds
              call xit("restart_write",-1)
            endif

            !--- Assign these values to a 2d array for the write to the restart
            vals2d(:,:) = reshape( vals(1:nwrds), (/nlon_ocn,nlat_ocn/) )
          endif

          !--- Write this field to the restart file
          ncstat =  nf90_inq_varid(rs_fid, trim(cpl_rs_var_name_ocn(idx1)), vid)
          call nc_error_handler(rs_fid, ncstat, 5)
          ncstat = nf90_put_var( rs_fid, vid, vals2d )
          call nc_error_handler(rs_fid, ncstat, 6, "writing "//trim(var_name))

          write(6,'(3a,i10,2a)')"restart_write: Wrote ",trim(cpl_rs_var_name_ocn(idx1)), &
                         " at time ",time_out," to ",trim(rs_file)
          call flush(6)
        enddo
        if ( associated(vals)   ) deallocate( vals )
        if ( associated(vals2d) ) deallocate( vals2d )
      endif

      !--- Close the restart file
      ncstat = nf90_close(rs_fid)
      call nc_error_handler(rs_fid, ncstat, 99)

    end subroutine restart_write

    !***************************************************************************
    !--- Assign initial values to a restart field, given its name
    !***************************************************************************
    subroutine restart_get_initial_values(vname, nx, ny, vals)
      character(*), intent(in)    :: vname
      integer(kind=4), intent(in) :: nx, ny
      real(kind=8), intent(inout) :: vals(nx,ny)

      !--- Local
      integer :: nwrds
      real(kind=8), pointer :: vals1d(:) => null()

      select case (trim(adjustl(vname)))
        case ("GT_atm")
          !--- Ground temperature
          vals(:,:) = 0.0

        case ("GC_atm")
          !--- Ground cover mask
          vals(:,:) = 0.0

        case ("SIC_atm")
          !--- Sea ice thickness
          vals(:,:) = 0.0

        case ("SICN_atm")
          !--- Sea ice fraction
          vals(:,:) = 0.0

        case ("SNO_atm")
          !--- Snow thickness
          vals(:,:) = 0.0

        case ("RES_atm")
          !--- Residual over lakes
          vals(:,:) = 0.0

        case ("OIceFrc")
          !--- Atmosphere-Ocean wind stress
          vals(:,:) = 0.0

        case ("O_OTaux1")
          !--- Atmosphere-Ocean wind stress
          vals(:,:) = 0.0

        case ("O_OTauy1")
          !--- Atmosphere-Ocean wind stress
          vals(:,:) = 0.0

        case ("O_OTaux2")
          !--- Atmosphere-Ocean wind stress
          !--- Set O_OTaux2 to O_OTaux1, if O_OTaux1 exists in storage
          !--- Otherwise initialize O_OTaux2 to zero
          call store_assign2d("O_OTaux1", ocn_gid, vals, nx, ny)

        case ("O_OTauy2")
          !--- Atmosphere-Ocean wind stress
          !--- Set O_OTauy2 to O_OTauy1, if O_OTauy1 exists in storage
          !--- Otherwise initialize O_OTauy2 to zero
          call store_assign2d("O_OTauy1", ocn_gid, vals, nx, ny)

        case ("O_ITaux1")
          !--- Atmosphere-Ice wind stress
          vals(:,:) = 0.0

        case ("O_ITauy1")
          !--- Atmosphere-Ice wind stress
          vals(:,:) = 0.0

        case ("O_ITaux2")
          !--- Atmosphere-Ice wind stress
          !--- Set O_ITaux2 to O_ITaux1, if O_ITaux1 exists in storage
          !--- Otherwise initialize O_ITaux2 to zero
          call store_assign2d("O_ITaux1", ocn_gid, vals, nx, ny)

        case ("O_ITauy2")
          !--- Atmosphere-Ice wind stress
          !--- Set O_ITauy2 to O_ITauy1, if O_ITauy1 exists in storage
          !--- Otherwise initialize O_ITauy2 to zero
          call store_assign2d("O_ITauy1", ocn_gid, vals, nx, ny)

        case ("O_QsrOce")
          !--- Solar heat flux over ocean
          vals(:,:) = 0.0

        case ("O_QsrIce")
          !--- Solar heat flux over ice
          vals(:,:) = 0.0

        case ("O_QsrMix")
          !--- Mixed solar heat flux over ocean and ice
          !--- Set O_QsrMix to O_QsrOce, if O_QsrOce exists in storage
          !--- Otherwise initialize O_QsrMix to zero
          call store_assign2d("O_QsrOce", ocn_gid, vals, nx, ny)

        case ("O_QnsOce")
          !--- Non-Solar heat flux over ocean
          vals(:,:) = 0.0

        case ("O_QnsIce")
          !--- Non-Solar heat flux over ice
          vals(:,:) = 0.0

        case ("O_QnsMix")
          !--- Mixed non-Solar heat flux over ocean and ice
          !--- Set O_QnsMix to O_QnsOce, if O_QnsOce exists in storage
          !--- Otherwise initialize O_QnsMix to zero
          call store_assign2d("O_QnsOce", ocn_gid, vals, nx, ny)

        case ("OTotRain")
          !--- Total liquid precipitation
          vals(:,:) = 0.0

        case ("OTotSnow")
          !--- Total solid precipitation
          vals(:,:) = 0.0

        case ("OTotEvap")
          !--- Total evaporation
          vals(:,:) = 0.0

        case ("OIceEvap")
          !--- Evaporation over ice
          vals(:,:) = 0.0

        case ("O_Wind10")
          !--- 10 meter wind
          vals(:,:) = 0.0

        case ("O_dQnsdT")
          !--- Non-solar sensitivity to temperature (dQns/dT)
          vals(:,:) = 0.0

        case ("O_Runoff")
          !--- Runoff
          vals(:,:) = 0.0

        case ("PMSL_atm")
          !--- mean sea level pressure (hPa)
          vals(:,:) = 1000.0

        case ("O_MSLP")
          !--- mean sea level pressure (atm)
          vals(:,:) = 1.0

        case ("CO2_atm")
          !--- CO2 concentration
          vals(:,:) = 288.0

        case ("O_AtmCO2")
          !--- CO2 concentration
          vals(:,:) = 288.0

        case ("CO2flx_atm")
          !--- CO2 flux
          vals(:,:) = 0.0

        case ("O_CO2FLX")
          !--- CO2 flux
          vals(:,:) = 0.0

        case ("OCalving")
          !--- Calving
          vals(:,:) = 0.0

        case ("O_TauMod")
          !--- Wind stress amplitude
          vals(:,:) = 0.0

        case ("OTopMlt")
          !--- top melt
          vals(:,:) = 0.0

        case ("OBotMlt")
          !--- bottom melt
          vals(:,:) = 0.0

        case default
          !--- Otherwise return a zero field
          vals(:,:) = 0.0
      end select

      !--- Clean up
      if ( associated(vals1d)   ) deallocate( vals1d )
      
    end subroutine restart_get_initial_values

    !***************************************************************************
    !--- Assign a 2D array with data from a storage variable
    !***************************************************************************
    subroutine store_assign2d(name, gid, vals2d, nx, ny)
      character(*), intent(in)    :: name
      integer, intent(in) :: gid
      real(kind=8), intent(inout) :: vals2d(:,:)
      integer(kind=4), intent(in) :: nx, ny

      !--- Local
      integer :: nwrds
      real(kind=8), pointer ::  vals1d(:) => null()

      nwrds = store_GetSize( trim(adjustl(name)), gid )
      if ( nwrds > 0 ) then
        !--- If nwrds > 0 then name exists in storage
        !--- Use name data to assign vals2d
        nwrds = store_readRecord(trim(adjustl(name)), gid, vals1d, idx=-1)
        if ( nwrds == nx*ny ) then
          vals2d(:,:) = reshape( vals1d(1:nwrds), (/nx,ny/) )
        else
          write(6,*)"store_assign2d: ",trim(name), &
              " contains no valid data. Initialize to zero."
          vals2d(:,:) = 0.0
        endif
      else
        write(6,*)"store_assign2d: ",trim(name), &
            " does not exist in storage. Initialize to zero."
        vals2d(:,:) = 0.0
      endif

      !--- Clean up
      if ( associated(vals1d) ) deallocate( vals1d )

    end subroutine store_assign2d

    !***************************************************************************
    !--- Create a CanCPL restart file
    !***************************************************************************
    subroutine restart_create(rs_file, rs_fid, time_string)
      character(*),    intent(in)  :: rs_file
      character(*),    intent(in)  :: time_string
      integer(kind=4), intent(out) :: rs_fid

      !--- Local
      integer (kind=4) :: ncstat, vid, idx
      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 :: idx1, grid_file_type
      character(32) :: var_name

      character(64) :: method = "restart_create"

      !--- 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(:,:) => null(), mid_lon_ocn(:,:) => null()
      real(kind=8), pointer ::  mid_lat_atm(:,:) => null(), mid_lat_ocn(:,:) => null()

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

      !--- Ensure the list of restart variables is defined
      call restart_init()

      !--- Create the output netcdf file in netcdf 4 format
      ncstat = nf90_create(rs_file, NF90_NETCDF4, rs_fid)
      call nc_error_handler(rs_fid, ncstat, 1, "Creating restart file "//trim(rs_file))

      !--- Write global attributes
      ncstat = nf90_put_att(rs_fid, NF90_GLOBAL, "restart_time", trim(time_string) )
      call nc_error_handler(rs_fid, ncstat, 2, "adding global attribute restart_time")

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

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

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

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

      if ( len_trim(cpl_runtime_repo_path) > 0 ) then
        ncstat = nf90_put_att(rs_fid, NF90_GLOBAL, "CanCPL_runtime_repo", trim(cpl_runtime_repo_path) )
        call nc_error_handler(rs_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(rs_fid, NF90_GLOBAL, "NEMO_build_commit", trim(nemo_commit) )
          call nc_error_handler(rs_fid, ncstat, 7, "adding global attribute NEMO_build_commit")
        endif

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

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

        if ( len_trim(nemo_build_cancpl_commit) > 0 ) then
          ncstat = nf90_put_att(rs_fid, NF90_GLOBAL, "NEMO_build_cancpl_commit", trim(nemo_build_cancpl_commit) )
          call nc_error_handler(rs_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(rs_fid, NF90_GLOBAL, "NEMO_build_cancpl_repo", trim(nemo_build_cancpl_repo_path) )
          call nc_error_handler(rs_fid, ncstat, 8, "adding global attribute NEMO_build_cancpl_repo")
        endif

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

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

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

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

      ncstat = nf90_def_dim(rs_fid, "nlat_atm", nlat_atm, lat_dim_atm_id)
      call nc_error_handler(rs_fid, ncstat, 7, "Defining lat dim atm")

      ncstat = nf90_def_dim(rs_fid, "nlon_ocn", nlon_ocn, lon_dim_ocn_id)
      call nc_error_handler(rs_fid, ncstat, 8, "Defining lon dim ocn")

      ncstat = nf90_def_dim(rs_fid, "nlat_ocn", nlat_ocn, lat_dim_ocn_id)
      call nc_error_handler(rs_fid, ncstat, 9, "Defining lat dim ocn")

      !--- Define ATM coordinate variables
      lon_atm_id = -1
      lat_atm_id = -1
      if ( len_trim(atm_grid_desc) > 0 ) then
        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)
          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)
          case default
            write(6,*)'restart_create: ',trim(atm_grid_desc), &
                '  Invalid grid_file_type = ',grid_file_type
            call xit("restart_create",-1)
        end select
        write(6,*)"restart_create: ATM grid_size, grid_dims = ",grid_size_atm,grid_dims_atm
        if ( associated(mid_lon_atm) ) then
          !--- Define the ATM lon coordinate as a 2d array
          ncstat = nf90_def_var(rs_fid, "lon_atm", NF90_DOUBLE, &
                                (/lon_dim_atm_id, lat_dim_atm_id/), lon_atm_id)
          call nc_error_handler(rs_fid, ncstat, 6, "defining lon_atm")
        endif
        if ( associated(mid_lat_atm) ) then
          !--- Define the ATM lat coordinate as a 2d array
          ncstat = nf90_def_var(rs_fid, "lat_atm", NF90_DOUBLE, &
                                (/lon_dim_atm_id, lat_dim_atm_id/), lat_atm_id)
          call nc_error_handler(rs_fid, ncstat, 7, "defining lat_atm")
        endif

      else
        !--- Extract atm lon from the internal atm grid definition
        call get_atm_lon(mid_lon_atm, overlap=.true.)

        !--- Define the ATM lon coordinate as a 2d array
        ncstat = nf90_def_var(rs_fid, "lon_atm", NF90_DOUBLE, &
                              (/lon_dim_atm_id, lat_dim_atm_id/), lon_atm_id)
        call nc_error_handler(rs_fid, ncstat, 6, "defining lon_atm")

        !--- Extract atm lat from the internal atm grid definition
        call get_atm_lat(mid_lat_atm)

        !--- Define the ATM lat coordinate as a 2d array
        ncstat = nf90_def_var(rs_fid, "lat_atm", NF90_DOUBLE, &
                              (/lon_dim_atm_id, lat_dim_atm_id/), lat_atm_id)
        call nc_error_handler(rs_fid, ncstat, 7, "defining lat_atm")

      endif

      !--- Define OCN coordinate variables
      lon_ocn_id = -1
      lat_ocn_id = -1
      if ( len_trim(ocn_grid_desc) > 0 ) then
        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)
          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)
          case default
            write(6,*)'restart_create: ',trim(ocn_grid_desc), &
                '  Invalid grid_file_type = ',grid_file_type
            call xit("restart_create",-1)
        end select
        write(6,*)"restart_create: OCN grid_size, grid_dims = ",grid_size_ocn,grid_dims_ocn
        if ( associated(mid_lon_ocn) ) then
          !--- Define the OCN lon coordinate as a 2d array
          ncstat = nf90_def_var(rs_fid, "lon_ocn", NF90_DOUBLE, &
                                (/lon_dim_ocn_id, lat_dim_ocn_id/), lon_ocn_id)
          call nc_error_handler(rs_fid, ncstat, 8, "defining lon_ocn")
        endif
        if ( associated(mid_lat_ocn) ) then
          !--- Define the OCN lat coordinate as a 2d array
          ncstat = nf90_def_var(rs_fid, "lat_ocn", NF90_DOUBLE, &
                                (/lon_dim_ocn_id, lat_dim_ocn_id/), lat_ocn_id)
          call nc_error_handler(rs_fid, ncstat, 9, "defining lat_ocn")
        endif

      else
        !--- Extract ocn lon from the internal ocn grid definition
        call get_ocn_lon(mid_lon_ocn, overlap=.true.)

        !--- Define the OCN lon coordinate as a 2d array
        ncstat = nf90_def_var(rs_fid, "lon_ocn", NF90_DOUBLE, &
                              (/lon_dim_ocn_id, lat_dim_ocn_id/), lon_ocn_id)
        call nc_error_handler(rs_fid, ncstat, 6, "defining lon_ocn")

        !--- Extract ocn lat from the internal ocn grid definition
        call get_ocn_lat(mid_lat_ocn)

        !--- Define the OCN lat coordinate as a 2d array
        ncstat = nf90_def_var(rs_fid, "lat_ocn", NF90_DOUBLE, &
                              (/lon_dim_ocn_id, lat_dim_ocn_id/), lat_ocn_id)
        call nc_error_handler(rs_fid, ncstat, 7, "defining lat_ocn")
      endif

      if ( nvar_atm > 0 ) then
        !--- Define AGCM restart variables
        do idx1=1,size(cpl_rs_var_name_atm)
          var_name = cpl_rs_var_name_atm(idx1)
          ncstat = nf90_def_var(rs_fid, trim(adjustl(var_name)), NF90_DOUBLE, &
                                (/lon_atm_id, lat_atm_id/), vid)
          call nc_error_handler(rs_fid, ncstat, 14, "defining "//trim(var_name))

          !--- Add variable specific meta data
          call nc_add_var_meta_data(var_name, rs_fid, vid)

          write(6,'(4a)')"restart_create: ",trim(cpl_rs_var_name_atm(idx1))," defined in ",trim(rs_file)
          call flush(6)
        enddo
      endif

      if ( nvar_ocn > 0 ) then
        !--- Define ocean restart variables
        do idx1=1,size(cpl_rs_var_name_ocn)
          var_name = cpl_rs_var_name_ocn(idx1)
          ncstat = nf90_def_var(rs_fid, trim(adjustl(var_name)), NF90_DOUBLE, &
                                (/lon_ocn_id, lat_ocn_id/), vid)
          call nc_error_handler(rs_fid, ncstat, 14, "defining "//trim(var_name))

          !--- Add variable specific meta data
          call nc_add_var_meta_data(var_name, rs_fid, vid)

          write(6,'(4a)')"restart_create: ",trim(cpl_rs_var_name_ocn(idx1))," defined in ",trim(rs_file)
          call flush(6)
        enddo
      endif

      !--- Get out of define mode
      ncstat = nf90_enddef(rs_fid)
      call nc_error_handler(rs_fid, ncstat, 16)

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

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

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

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

    end subroutine restart_create

    !***************************************************************************
    !--- Read certain variables from a NEMO restart file
    !***************************************************************************
    subroutine read_nemo_restart(name, gid, found)
      character(*), intent(in) :: name
      integer, intent(in) :: gid
      logical, intent(out) :: found

      !--- Local
      character(128) :: nemo_ocn_rs, nemo_ice_rs
      integer :: iu
      logical :: exists
      integer(4) :: fid, vid, ncstat, nwrds, nx, nxid, ny, nyid, nz, nzid
      real(8), pointer :: src_data(:) => null(), dst_data(:) => null()
      real(8), pointer :: awrk2d(:,:) => null(), bwrk2d(:,:) => null()
      real(8), pointer :: wrk3d(:,:,:) => null()

      !--- These are the constant values for sea ice and snow density
      !--- that are used by lim2
      real(8), parameter :: rhoic=900.0_8  !--- Kg m-3
      real(8), parameter :: rhosn=330.0_8  !--- Kg m-3

      !--- The name of the NEMO ocean and ice restart files is currently
      !--- hard coded in nemo_prelude to be the following.
      !--- This may change.
      nemo_ocn_rs = "NEMO_ocn_global_rsin.nc"
      nemo_ice_rs = "NEMO_ice_global_rsin.nc"

      found = .false.
      select case( trim(adjustl(name)) )
        case("O_SSTSST", "GT_atm")
          !--- Sea surface temperature (K)
          !--- This is found in the NEMO ocean restart file
          inquire(file=trim(nemo_ocn_rs), exist=exists)
          if ( exists ) then

            !--- Read the 3D ocean temperature from the ocean restart file
            !--- NOTE: tn contains temperature in Celcius
            call ncdf_open(fid, nemo_ocn_rs, "r")
            call ncdf_read_var(fid, "tn", wrk3d, nt=1)
            call ncdf_close(fid)

            if ( nlon_ocn /= size(wrk3d,dim=1) .or. nlat_ocn /= size(wrk3d,dim=2) ) then
              write(6,*)"read_nemo_restart: Dimension size mismatch for ",trim(name)
              write(6,*)"                   Reading tn from file ",trim(nemo_ocn_rs)
              write(6,*)"      read nx,ny: ",size(wrk3d,dim=1),size(wrk3d,dim=2)
              write(6,*)"  expected nx,ny: ",nlon_ocn,nlat_ocn
              call xit("read_nemo_restart",-1)
            endif

            if ( associated(src_data) ) deallocate(src_data)
            allocate( src_data(nlon_ocn*nlat_ocn) )
            !--- Level 1 of tn is SST
            src_data = reshape(wrk3d(1:nlon_ocn,1:nlat_ocn,1), (/nlon_ocn*nlat_ocn/))

            !--- Convert from Celcius to Kelvin
            src_data = src_data + 273.15_8

            !--- Add this field to internal storage
            if ( trim(adjustl(name)) .eq. "GT_atm" ) then
              !--- Interpolate tn, read from the NEMO restart file, onto the agcm grid
              call regrid_field_cv("O_SSTSST", src_data, "GT_atm", dst_data)

#ifdef DBG_bfb
strng = " "
write(strng,'(2a)')"restart_var_interp_",trim(name)
call dump_farray(trim(strng),dst_data,idx=cpl_elapsed_time_secs)
#endif

              !--- Store the interpolated field
              call store_AppendRecord("GT_atm", gid, dst_data, create=.true.)
            else
              !--- Store tn, as it was read from the NEMO restart file
              call store_AppendRecord("O_SSTSST", gid, src_data, create=.true.)
            endif

            found = .true.
          else
            write(6,*)"read_nemo_restart: Missing ocean restart file ",trim(nemo_ocn_rs)
          endif

        case("O_TepIce")
          !--- Surface temperature of sea ice
          !--- This is found in the NEMO ice restart file
          inquire(file=trim(nemo_ice_rs), exist=exists)
          if ( exists ) then

            !--- Read the 2D ice surface temperature from the ice restart file
            call ncdf_open(fid, nemo_ice_rs, "r")
            call ncdf_read_var(fid, "sist", awrk2d, nt=1)
            call ncdf_close(fid)

            if ( nlon_ocn /= size(awrk2d,dim=1) .or. nlat_ocn /= size(awrk2d,dim=2) ) then
              write(6,*)"read_nemo_restart: Dimension size mismatch for ",trim(name)
              write(6,*)"                   Reading sist from file ",trim(nemo_ice_rs)
              write(6,*)"      read nx,ny: ",size(awrk2d,dim=1),size(awrk2d,dim=2)
              write(6,*)"  expected nx,ny: ",nlon_ocn,nlat_ocn
              call xit("read_nemo_restart",-1)
            endif

            if ( associated(src_data) ) deallocate(src_data)
            allocate( src_data(nlon_ocn*nlat_ocn) )
            !--- Assign sist to a 1D array
            src_data = reshape(awrk2d(1:nlon_ocn,1:nlat_ocn), (/nlon_ocn*nlat_ocn/))

            !--- Add this field to internal storage
            !--- Store sist, as it was read from the NEMO restart file
            call store_AppendRecord("O_TepIce", gid, src_data, create=.true.)

            found = .true.
          else
            write(6,*)"read_nemo_restart: Missing ice restart file ",trim(nemo_ice_rs)
          endif

        case("OIceFrc", "SICN_atm")
          !--- Sea ice fraction
          !--- This is found in the NEMO ice restart file
          inquire(file=trim(nemo_ice_rs), exist=exists)
          if ( exists ) then

            !--- Read the 2D lead fraction from the ice restart file
            call ncdf_open(fid, nemo_ice_rs, "r")
            call ncdf_read_var(fid, "frld", awrk2d, nt=1)
            call ncdf_close(fid)

            if ( nlon_ocn /= size(awrk2d,dim=1) .or. nlat_ocn /= size(awrk2d,dim=2) ) then
              write(6,*)"read_nemo_restart: Dimension size mismatch for ",trim(name)
              write(6,*)"                   Reading frld from file ",trim(nemo_ice_rs)
              write(6,*)"      read nx,ny: ",size(awrk2d,dim=1),size(awrk2d,dim=2)
              write(6,*)"  expected nx,ny: ",nlon_ocn,nlat_ocn
              call xit("read_nemo_restart",-1)
            endif

            if ( associated(src_data) ) deallocate(src_data)
            allocate( src_data(nlon_ocn*nlat_ocn) )
            !--- Assign frld to a 1D array
            src_data = reshape(awrk2d(1:nlon_ocn,1:nlat_ocn), (/nlon_ocn*nlat_ocn/))

            !--- Convert frld to ice fraction
            src_data = 1.0_8 - src_data

            !--- Add this field to internal storage
            if ( trim(adjustl(name)) .eq. "SICN_atm" ) then
              !--- Interpolate 1-frld, read from the NEMO restart file, onto the agcm grid
              call regrid_field_cv("OIceFrc", src_data, "SICN_atm", dst_data)

#ifdef DBG_bfb
strng = " "
write(strng,'(2a)')"restart_var_interp_",trim(name)
call dump_farray(trim(strng),dst_data,idx=cpl_elapsed_time_secs)
#endif

              !--- The interpolation can result in values outside the range 0-1
              where (dst_data<0.0_8) dst_data = 0.0_8
              where (dst_data>1.0_8) dst_data = 1.0_8

#ifdef DBG_bfb
strng = " "
write(strng,'(2a)')"restart_var_interp2_",trim(name)
call dump_farray(trim(strng),dst_data,idx=cpl_elapsed_time_secs)
#endif

              !--- Store the interpolated field
              call store_AppendRecord("SICN_atm", gid, dst_data, create=.true.)
            else
              !--- Store 1-frld, as it was read from the NEMO restart file
              call store_AppendRecord("OIceFrc", gid, src_data, create=.true.)
            endif

            found = .true.
          else
            write(6,*)"read_nemo_restart: Missing ice restart file ",trim(nemo_ice_rs)
          endif

        case("OIceTck", "SIC_atm")
          !--- Sea Ice Water Equivalent (kg m-2)
          !--- This is found in the NEMO ice restart file
          inquire(file=trim(nemo_ice_rs), exist=exists)
          if ( exists ) then

            !--- Read the 2D ice thickness (hicif) and lead fraction (frld)
            !--- from the ice restart file
            call ncdf_open(fid, nemo_ice_rs, "r")
            call ncdf_read_var(fid, "hicif", awrk2d, nt=1)
            call ncdf_read_var(fid,  "frld", bwrk2d, nt=1)
            call ncdf_close(fid)

            if ( nlon_ocn /= size(awrk2d,dim=1) .or. nlat_ocn /= size(awrk2d,dim=2) ) then
              write(6,*)"read_nemo_restart: Dimension size mismatch for ",trim(name)
              write(6,*)"                   Reading hicif from file ",trim(nemo_ice_rs)
              write(6,*)"      read nx,ny: ",size(awrk2d,dim=1),size(awrk2d,dim=2)
              write(6,*)"  expected nx,ny: ",nlon_ocn,nlat_ocn
              call xit("read_nemo_restart",-1)
            endif

            if ( nlon_ocn /= size(bwrk2d,dim=1) .or. nlat_ocn /= size(bwrk2d,dim=2) ) then
              write(6,*)"read_nemo_restart: Dimension size mismatch for ",trim(name)
              write(6,*)"                   Reading frld from file ",trim(nemo_ice_rs)
              write(6,*)"      read nx,ny: ",size(bwrk2d,dim=1),size(bwrk2d,dim=2)
              write(6,*)"  expected nx,ny: ",nlon_ocn,nlat_ocn
              call xit("read_nemo_restart",-1)
            endif

            !--- convert sea ice thickness, hicif (m), to sea ice water equivalent (kg m-2)
            !--- rhoic is the constant ice density that is used in LIM2
            awrk2d = rhoic * awrk2d * (1.0_8-bwrk2d)

            if ( associated(src_data) ) deallocate(src_data)
            allocate( src_data(nlon_ocn*nlat_ocn) )
            !--- Assign frld to a 1D array
            src_data = reshape(awrk2d(1:nlon_ocn,1:nlat_ocn), (/nlon_ocn*nlat_ocn/))

            !--- Add this field to internal storage
            if ( trim(adjustl(name)) .eq. "SIC_atm" ) then
              !--- Interpolate 1-frld, read from the NEMO restart file, onto the agcm grid
              call regrid_field_cv("OIceTck", src_data, "SIC_atm", dst_data)

#ifdef DBG_bfb
strng = " "
write(strng,'(2a)')"restart_var_interp_",trim(name)
call dump_farray(trim(strng),dst_data,idx=cpl_elapsed_time_secs)
#endif

              !--- Store the interpolated field
              call store_AppendRecord("SIC_atm", gid, dst_data, create=.true.)
            else
              !--- Store 1-frld, as it was read from the NEMO restart file
              call store_AppendRecord("OIceTck", gid, src_data, create=.true.)
            endif

            found = .true.
          else
            write(6,*)"read_nemo_restart: Missing ice restart file ",trim(nemo_ice_rs)
          endif

        case("OSnwTck", "SNO_atm")
          !--- Snow Water Equivalent over Sea Ice (kg m-2)
          !--- This is found in the NEMO ice restart file
          inquire(file=trim(nemo_ice_rs), exist=exists)
          if ( exists ) then

            !--- Read the 2D snow thickness (hsnif) and lead fraction (frld)
            !--- from the ice restart file
            call ncdf_open(fid, nemo_ice_rs, "r")
            call ncdf_read_var(fid, "hsnif", awrk2d, nt=1)
            call ncdf_read_var(fid,  "frld", bwrk2d, nt=1)
            call ncdf_close(fid)

            if ( nlon_ocn /= size(awrk2d,dim=1) .or. nlat_ocn /= size(awrk2d,dim=2) ) then
              write(6,*)"read_nemo_restart: Dimension size mismatch for ",trim(name)
              write(6,*)"                   Reading hsnif from file ",trim(nemo_ice_rs)
              write(6,*)"      read nx,ny: ",size(awrk2d,dim=1),size(awrk2d,dim=2)
              write(6,*)"  expected nx,ny: ",nlon_ocn,nlat_ocn
              call xit("read_nemo_restart",-1)
            endif

            if ( nlon_ocn /= size(bwrk2d,dim=1) .or. nlat_ocn /= size(bwrk2d,dim=2) ) then
              write(6,*)"read_nemo_restart: Dimension size mismatch for ",trim(name)
              write(6,*)"                   Reading frld from file ",trim(nemo_ice_rs)
              write(6,*)"      read nx,ny: ",size(bwrk2d,dim=1),size(bwrk2d,dim=2)
              write(6,*)"  expected nx,ny: ",nlon_ocn,nlat_ocn
              call xit("read_nemo_restart",-1)
            endif

            !--- convert snow thickness, hsnif (m), to snow water equivalent (kg m-2)
            !--- rhosn is the constant snow density that is used in LIM2
            awrk2d = rhosn * awrk2d * (1.0_8-bwrk2d)

            if ( associated(src_data) ) deallocate(src_data)
            allocate( src_data(nlon_ocn*nlat_ocn) )
            !--- Assign frld to a 1D array
            src_data = reshape(awrk2d(1:nlon_ocn,1:nlat_ocn), (/nlon_ocn*nlat_ocn/))

            !--- Add this field to internal storage
            if ( trim(adjustl(name)) .eq. "SNO_atm" ) then
              !--- Interpolate data read from the NEMO restart file, onto the agcm grid
              call regrid_field_cv("OSnwTck", src_data, "SNO_atm", dst_data)

#ifdef DBG_bfb
strng = " "
write(strng,'(2a)')"restart_var_interp_",trim(name)
call dump_farray(trim(strng),dst_data,idx=cpl_elapsed_time_secs)
#endif

              !--- Store the interpolated field
              call store_AppendRecord("SNO_atm", gid, dst_data, create=.true.)
            else
              !--- Store data read from the NEMO restart file on the NEMO grid
              call store_AppendRecord("OSnwTck", gid, src_data, create=.true.)
            endif

            found = .true.
          else
            write(6,*)"read_nemo_restart: Missing ice restart file ",trim(nemo_ice_rs)
          endif

      end select

#ifdef DBG_bfb
if ( found) then
  strng = " "
  write(strng,'(2a)')"restart_var_nemo_",trim(name)
  call dump_farray(trim(strng),src_data,idx=cpl_elapsed_time_secs)
endif
#endif

      !--- Clean up
      if ( associated(src_data) ) deallocate(src_data)
      if ( associated(dst_data) ) deallocate(dst_data)
      if ( associated(awrk2d) ) deallocate(awrk2d)
      if ( associated(bwrk2d) ) deallocate(bwrk2d)
      if ( associated(wrk3d) ) deallocate(wrk3d)

    end subroutine read_nemo_restart

end module restart

