module atm_f_data

  use com_cpl
  use ccc_io
  use store
  use ncdf
  use strings

  implicit none

  private

  public :: read_atm_forcing_data

  !--- Day of the year corresponding to the start of each month
  integer, parameter, dimension(12) :: som = &
    (/ 1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335 /)

  !--- Day of the year corresponding to the end of each month
  integer, parameter, dimension(12) :: eom = &
    (/31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 /)

  !--- These are ibuf2 values read from monthly climatology files
  !--- corresponding to mid month days
  integer, parameter, dimension(12) :: mmd = &
    (/16,46,75,106,136,167,197,228,259,289,320,350/)

  !--- Keep track of the first year/day for which boundary conditions are read
  integer, save :: specified_bc_start_year = -1
  integer, save :: specified_bc_start_day = -1

  !--- The names of netcdf file containing AGCM boundary condition data (GT,SIC,SICN)
  character(len=256), save :: gt_file_name, sic_file_name, sicn_file_name, gtano_file_name
  character(len=256), save :: clim_file_name

  !--- Used in calls to netcdf access routines
  integer(kind=4), save :: fid_gt, fid_sic, fid_sicn, fid_clim, fid_gtano, ncvid, ncstat
  integer(kind=4), save :: vid_gt, vid_sic, vid_sicn, vid_clim, vid_gtano
  real(kind=8), save, pointer :: vals_2dr8(:,:)
  real(kind=8), save :: gsize

  !--- Store the values of the time coordinate (record dimension) found in the
  !--- AGCM boundary condition netcdf data files
  real(kind=8), save, pointer :: time_vals(:)
  integer, save :: ntime = 0
  character(len=64), save :: time_units
  integer, save :: time_origin_ymd_hms(6)
  integer, save :: time_idx
  !--- Store the values of the time coordinate (record dimension) found in GTano file
  real(kind=8), save, pointer :: time_vals_gtano(:)
  integer, save :: ntime_gtano = 0
  character(len=64), save :: time_units_gtano
  integer, save :: time_origin_ymd_hms_gtano(6)
  integer, save :: time_idx_gtano

  !--- These 2 logical variables control whether or not the AGCM boundary
  !--- condition data is in CCCma binary format or netcdf.
  !--- They can both be true in which case the values read in at run time are
  !--- compared to find the min, max and relative difference but the netcdf data
  !--- are the values that are actually used.
  !--- It is anticipated that netcdf will be the format that is used exclusively.
  !--- The ability to read CCCma binary data was only used to verify that both
  !--- formats are bit identical when the netcdf file is converted from the
  !--- CCCma binary data prior to execution.
  logical, save :: atm_bc_from_netcdf = .true.
  logical, save :: atm_bc_from_cccma = .false.
  logical, save :: gtano_from_netcdf = .false.

  contains

    !***************************************************************************
    !--- Given the year and month return the index in the time_vals array that
    !--- represents that year and mid month day
    !***************************************************************************
    function time_slice(year, month) result(indx)
      integer :: indx
      integer, intent(in) :: year, month

      !--- Local
      real(kind=8) :: days
      integer :: idx

      indx = -1

      !--- Determine the number of days since time origin that represents the
      !--- current year and month
      if ( specified_bc == 1 .or. specified_bc == 2 ) then
        !--- This is climatological data
        !--- This data has been converted from a CCCma AGCM "AN" file which
        !--- used an ibuf2 format which was mid month day of the year
        !--- Therefore use mid month day of year
        days = mmd(month)
      else
        !--- Otherwise assume multi year data
        !--- This data has been converted from CCCma binary files that used an
        !--- ibuf2 format of YYYYMMDDHH in which DDHH was always 0100
        !--- Because of this, the day associated with any month is always the
        !--- first day of that month
        days = (year - time_origin_ymd_hms(1))*365 + som(month) - 1
      endif

      do idx=1,size(time_vals)
        if ( abs(time_vals(idx)-days) < 0.01 ) then
          indx = idx
          exit
        endif
      enddo

      if ( indx < 0 ) then
        write(6,*)"time_slice: Unable to determine time index for year, mon = ", &
                  year,month
        call xit("time_slice",-1)
      endif

    end function time_slice

    function time_slice_gtano(year, month) result(indx)
      integer :: indx
      integer, intent(in) :: year, month

      !--- Local
      real(kind=8) :: days
      integer :: idx

      indx = -1

      !--- Determine the number of days since time origin that represents the
      !--- current year and month

      !--- This data has been converted from CCCma binary files that used an
      !--- ibuf2 format of YYYYMMDDHH in which DDHH was always 0100
      !--- Because of this, the day associated with any month is always the
      !--- first day of that month
      days = (year - time_origin_ymd_hms_gtano(1))*365 + som(month) - 1

      do idx=1,size(time_vals_gtano)
        if ( abs(time_vals_gtano(idx)-days) < 0.01 ) then
          indx = idx
          exit
        endif
      enddo

      if ( indx < 0 ) then
        write(6,*)"time_slice_gtano: Unable to determine time index for year, mon = ", &
                  year,month
        call xit("time_slice_gtano",-1)
      endif

    end function time_slice_gtano

    !***************************************************************************
    !--- Read atm forcing data from a file
    !***************************************************************************
    subroutine read_atm_forcing_data()

      !--- Local
      integer :: nwrds, gid, verbose=6
      real(kind=8), pointer :: wrk1(:)
      integer(kind=8) :: store_time

      !--- first and second dimension of input arrays
      integer :: lon, lat

      !--- Fields to be returned from read_spec_bc
      real, allocatable :: gt(:,:), sic(:,:), sicn(:,:)

      !--- Field to be input, possibly modified, and returned from read_spec_bc
      real, allocatable :: sno(:,:), flnd(:,:)

      !--- current year and day of year at which to interpolate
      integer :: year, day

      !--- Toggle to determine if we abort when fields are not found in storage
      logical :: abort_if_missing=.false.

      logical :: exists, found

      nullify(wrk1)

      !--- Validate the integer flag indicating the type of boundary data
      if ( specified_bc < 1 .or. specified_bc > 4 ) then
        write(6,*)"read_atm_forcing_data: specified_bc=",specified_bc," is out of range."
        call xit ("read_atm_forcing_data",-1)
      endif

      !--- Ensure the file name is supplied for the forcing data file
      if (len_trim(specified_bc_file).le.0) then
        write(6,*)"read_atm_forcing_data: specified_bc_file is blank."
        call xit("read_atm_forcing_data",-2)
      endif

      if ( trim(adjustl(specified_bc_file)) .eq. "AN" .or. &
           trim(adjustl(specified_bc_file)) .eq. "AN.nc" ) then
        !--- The "AN" file only contains monthly averaged climatologies
        if (specified_bc .ne. 1) then
          write(6,*)"read_atm_forcing_data: Attempting to read AN file with specified_bc=", &
              specified_bc
          call xit ("read_atm_forcing_data",-3)
        endif
      endif

      !--- Define netcdf file names based on the prefix specified_bc_file
      gt_file_name   = trim(specified_bc_file)//"_gt.nc"
      sic_file_name  = trim(specified_bc_file)//"_sic.nc"
      sicn_file_name = trim(specified_bc_file)//"_sicn.nc"
      clim_file_name = trim(specified_bc_file)//".nc"
      gtano_file_name = trim(specified_bc_file)//"_gtano.nc"

      if ( atm_bc_from_netcdf ) then

        !--- Check if the GTano netcdf data file exists in the current working dir
        inquire(file=trim(gtano_file_name),exist=exists)
        if ( exists ) then
          write(6,*)"read_atm_forcing_data: GTANO file name is found",trim(gtano_file_name)
          gtano_from_netcdf = .true.
          !--- Determine the variable ID for GTano in this file
          vid_gtano = ncdf_get_varid(gtano_file_name, "GT_atm")
          call ncdf_get_time(gtano_file_name, time_vals_gtano, ntime=ntime_gtano, units=time_units_gtano)
          write(6,*)"read_atm_forcing_data: Time units found in ",trim(gtano_file_name)," are ",trim(time_units_gtano)
          !--- Break these units down into year, month, day, hour, minute, second and store
          !--- these values in the global array time_origin_ymd_hms
          call parse_time_units(time_units_gtano)
          time_origin_ymd_hms_gtano = time_origin_ymd_hms
        else
          write(6,*)"read_atm_forcing_data: GTANO file name is not found",trim(gtano_file_name)
          gtano_from_netcdf = .false.
        endif

        if ( specified_bc == 1 .or. specified_bc == 2 ) then
          !--- Climatological data

          !--- Ensure that the netcdf bc data files exist in the current working dir
          if ( .not. locate_file(clim_file_name) ) then
            write(6,*)"read_atm_forcing_data: ",trim(clim_file_name)," is missing."
            call xit("read_atm_forcing_data",-4)
          endif

          !--- Determine the variable ID for GT_atm, SIC_atm, SICN_atm in this file
          vid_gt   = ncdf_get_varid(clim_file_name, "GT_atm")
          vid_sic  = ncdf_get_varid(clim_file_name, "SIC_atm")
          vid_sicn = ncdf_get_varid(clim_file_name, "SICN_atm")

          !--- Define an array containing time values found in the netcdf data file
          if ( .not. associated(time_vals) ) then
            call ncdf_get_time(clim_file_name, time_vals, ntime=ntime, units=time_units)
            if ( ntime == 0 ) then
              write(6,*)"read_spec_bc: No time slices are present in ",trim(clim_file_name)
              call xit("read_atm_forcing_data",-32)
            endif
            write(6,*)"read_atm_forcing_data: Time units found in ",trim(clim_file_name)," are ",trim(time_units)
            !--- Break these units down into year, month, day, hour, minute, second and store
            !--- these values in the global array time_origin_ymd_hms
            call parse_time_units(time_units)
          endif
        endif

        if ( specified_bc == 3 .or. specified_bc == 4 ) then
          !--- Multi year data

          !--- Ensure that the GT netcdf bc data file exists in the current working dir
          if ( .not. locate_file(gt_file_name) ) then
            write(6,*)"read_atm_forcing_data: ",trim(gt_file_name)," is missing."
            call xit("read_atm_forcing_data",-5)
          endif

          !--- Determine the variable ID for GT_atm in this file
          vid_gt = ncdf_get_varid(gt_file_name, "GT_atm")

          !--- Ensure that the SIC netcdf bc data file exists in the current working dir
          if ( .not. locate_file(sic_file_name) ) then
            write(6,*)"read_atm_forcing_data: ",trim(sic_file_name)," is missing."
            call xit("read_atm_forcing_data",-6)
          endif

          !--- Determine the variable ID for SIC_atm in this file
          vid_sic = ncdf_get_varid(sic_file_name, "SIC_atm")

          !--- Ensure that the SICN netcdf bc data file exists in the current working dir
          if ( .not. locate_file(sicn_file_name) ) then
            write(6,*)"read_atm_forcing_data: ",trim(sicn_file_name)," is missing."
            call xit("read_atm_forcing_data",-7)
          endif

          !--- Determine the variable ID for SICN_atm in this file
          vid_sicn = ncdf_get_varid(sicn_file_name, "SICN_atm")

          !--- Define an array containing time values found in the netcdf data files
          !--- Assume all three files contain the same time slices
          if ( .not. associated(time_vals) ) then
            call ncdf_get_time(gt_file_name, time_vals, ntime=ntime, units=time_units)
            if ( ntime == 0 ) then
              write(6,*)"read_spec_bc: No time slices are present in ",trim(gt_file_name)
              call xit("read_atm_forcing_data",-33)
            endif
            write(6,*)"read_atm_forcing_data: Time units found in ",trim(gt_file_name)," are ",trim(time_units)
            !--- Break these units down into year, month, day, hour, minute, second and store
            !--- these values in the global array time_origin_ymd_hms
            call parse_time_units(time_units)
          endif
        endif
      endif

      lon = nlon_a + olap_a
      lat = nlat_a

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

      !--- Determine the current year and day of the year
      !--- read_spec_bc will only interpolate down to a day
      year = cpl_year
      day  = cpl_day_of_year

      if ( specified_bc_start_year < 0 ) specified_bc_start_year = year
      if ( specified_bc_start_day < 0 )  specified_bc_start_day  = day

      !--- Allocate arrays passed to read_spec_bc
      allocate( gt(lon,lat), sic(lon,lat), sicn(lon,lat), &
               sno(lon,lat), flnd(lon,lat) )

      !--- Initialize arrays passed to read_spec_bc
      gt   = 0.0
      sic  = 0.0
      sicn = 0.0

      !--- Initialize SNO
      nwrds = -1
      if ( store_GetIndex("SNO_atm",anygid) > 0 ) then
        !--- Get SNO from internal storage
        !--- This will be the value sent from the AGCM on the previous time step
        !--- NOTE store_readRecord will return nwrds=-1 if storage for SNO_atm is
        !---      initialized but does not contain any records
        nwrds = store_readRecord("SNO_atm", anygid, wrk1, idx=-1)
      endif

      if ( nwrds < 0 ) then
        !--- The field was not found in storage
        if ( abort_if_missing ) then
          write(6,*)"read_atm_forcing_data: SNO_atm is missing from storage."
          call xit("read_atm_forcing_data",-4)
        else
          nwrds = lon*lat
          if ( associated(wrk1) ) deallocate(wrk1)
          allocate( wrk1(nwrds) )
          wrk1=0.0
        endif
      endif

      if ( nwrds == lon*lat ) then
        sno(1:lon,1:lat) = reshape( wrk1, (/lon,lat/) )
      else
        write(6,*)"read_atm_forcing_data: Size of SNO_atm is incorrect."
        write(6,*)"       expected size of SNO_atm = ",lon*lat
        write(6,*)"          found size of SNO_atm = ",nwrds
        call xit("read_atm_forcing_data",-5)
      endif

      !--- Initialize FLND from atm_fland
      if ( associated(atm_fland) ) then
        if ( size(atm_fland) == lon*lat ) then
          flnd(1:lon,1:lat) = reshape( atm_fland, (/lon,lat/) )
        else
          write(6,*)"read_atm_forcing_data: size mismatch for atm_fland"
          write(6,*)"Expected size: ",lon*lat
          write(6,*)"  Actual size: ",size(atm_fland)
          call xit("read_atm_forcing_data",-7)
        endif
      else
        write(6,*)"read_atm_forcing_data: atm_fland is not allocated."
        call xit("read_atm_forcing_data",-8)
      endif

      !--- Read the specified boundary conditions from a file and
      !--- interpolate to the current day
      call read_spec_bc(gt,sic,sicn,sno,flnd,lon,lat,year,day, &
                        specified_bc,specified_bc_file,verbose)

      !--- Write the fields returned from read_spec_bc to internal storage
      !--- These fields are associated with the ocn group because if the ocean was
      !--- coupled they would have been interpolated from data that was received
      !--- from the ocean
      if ( associated(wrk1) ) deallocate(wrk1)
      allocate( wrk1(lon*lat) )

      !--- Ground temperature
      wrk1(1:lon*lat) = reshape( gt, (/lon*lat/) )
      gid = filter_group_ID("GT_atm", atm_gid)
      call store_AppendRecord("GT_atm", gid, wrk1, create=.true.)

      !--- Sea ice fraction
      wrk1(1:lon*lat) = reshape( sicn, (/lon*lat/) )
      gid = filter_group_ID("SICN_atm", atm_gid)
      call store_AppendRecord("SICN_atm", gid, wrk1, create=.true.)

      !--- Sea ice thickness over sea ice (aka sea ice water equivalent)
      wrk1(1:lon*lat) = reshape( sic, (/lon*lat/) )
      gid = filter_group_ID("SIC_atm", atm_gid)
      call store_AppendRecord("SIC_atm", gid, wrk1, create=.true.)

      !--- Snow thickness over sea ice (aka snow water equivalent)
      wrk1(1:lon*lat) = reshape( sno, (/lon*lat/) )
      gid = filter_group_ID("SNO_atm", atm_gid)
      call store_AppendRecord("SNO_atm", gid, wrk1, create=.true.)

      if (verbose.gt.10) then
        write(6,*)"read_atm_forcing_data: Fields to be saved on internal storage"
        call store_ListStorage()
      endif

      !--- Clean up
      deallocate( gt, sic, sicn, sno, flnd )
      if ( associated(wrk1) ) deallocate(wrk1)

    end subroutine read_atm_forcing_data

    !***********************************************************************
    !--- Read specified boundary conditions for GT,SIC and SICN from a file.
    !--- This subroutine is from the file cpl_read_spec_bc.F which was/is
    !--- part of the CanESM4.2 ocean code
    !--- Larry Solheim ...Jun,2016
    !***********************************************************************
    subroutine read_spec_bc(gt,sic,sicn,sno,flnd,lon,lat,year,day, &
                            specified_bc,specified_bc_file,verbose)
      !--- Input/Output

      !--- Fields to be returned
      real :: gt(lon,lat), sic(lon,lat), sicn(lon,lat)

      !--- Field to be input, possibly modified, and returned
      real :: sno(lon,lat), flnd(lon,lat)

      !--- first and second dimension of input arrays
      integer :: lon, lat

      !--- current year and day of year at which to interpolate
      integer :: year, day

      !--- An integer flag to indicate the type of boundary data
      !--- to be read from the file
      !--- specified_bc = 1 ...Monthly climatology
      !--- specified_bc = 2 ...Daily climatology
      !--- specified_bc = 3 ...multi-year monthly data
      !--- specified_bc = 4 ...multi-year daily data
      integer :: specified_bc

      !--- Name of a file containing the boundary data
      character*(*) :: specified_bc_file

      !--- provide debugging/diagnostic info when verbose > 0
      integer :: verbose

      !--- local
      logical :: remote, ok, check_ice
      character(len=512) :: strng
      character(len=512) :: caller_file_name = __FILE__
      integer :: caller_line
      integer :: year1, year2, mon1, mon2
      integer :: i, j, k, month, dom, year_offset
      integer :: nx1, ny1, nx2, ny2
      integer :: nf21=21, nf22=22, nf23=23
      integer :: useyear, usemonth, timestmp
      real :: clim_time(3,1:14), tmp_time(3), gavg, sicmax
      real :: w1, w2, den, dd
      integer :: nc4to8
      external   nc4to8

      integer, parameter :: maxx=131328  !--- 256x513
      integer :: ibuf(8),idat(maxx)
      common /icom/ ibuf,idat

      !--- Extra variables required for calls to PHSIL8
      real :: gtrol(lon,lat), sicrol(lon,lat), sicnrol(lon,lat)
      real :: gtrom(lon,lat), sicrom(lon,lat), sicnrom(lon,lat)
      real :: gtano(lon,lat), gtanom(lon,lat), gtanol(lon,lat)
      integer :: ilg, il1, il2, ib2rom, ib2rol, iyear, iday, mday
      integer :: curr_mon
      real :: delt, gmt

      !--- SICN_CRT is the minimum value of SICN for which a grid box
      !--- is considered to be ice covered.
      !--- Historically, in the agcm SICN_CRT=0.175 for climatologies whereas
      !--- for multi-year boundary conditions SICN_CRT=0.15
      !--- In coupled runs, always fix SICN_CRT=0.15
      real, parameter :: sicn_crt = 0.15

      !--- Dummy arguments passed back from SPWCON11
      real :: fvort, pii

      !--- Common block used to define GTFSW in current routine
      !--- These values are defined by the call to SPWCON11
      real ::         csno, cpack, gtfsw, rkhi, sbc, snomax
      common /param3/ csno, cpack, gtfsw, rkhi, sbc, snomax

      !.......Check on input values
      if (year.lt.1) then
        write(6,*)'READ_SPEC_BC: year=',year,' is out of range.'
        call xit("READ_SPEC_BC",-1)
      endif
      if (day.lt.1 .or. day.gt.365) then
        write(6,*)'READ_SPEC_BC: day=',day,' is out of range.'
        call xit("READ_SPEC_BC",-2)
      endif
      if (specified_bc.lt.1 .or. specified_bc.gt.4) then
        write(6,*)'READ_SPEC_BC: specified_bc=',specified_bc &
          ,' is out of range.'
         call xit("READ_SPEC_BC",-3)
      endif

      !--- Initialize common blocks used by PHSIL8
      call spwcon11 (fvort,pii)

      if (verbose.gt.10) then
        strng=' '
        write(strng,'(a,6g16.8)') &
          ' csno, cpack, gtfsw, rkhi, sbc, snomax: ', &
            csno, cpack, gtfsw, rkhi, sbc, snomax
        caller_line      = __LINE__
        call dbg_out(strng, caller_file_name, caller_line)
      endif

      !.......................................................................
      !.......Read GT/SIC/SICN from different files, depending on specified_bc
      !.......................................................................

      if (specified_bc .eq. 1) then
        !-------------------------------------------------------
        !--- Monthly averaged climatological boundary conditions
        !-------------------------------------------------------

        if (verbose.gt.0) then
          write(strng,'(2a)') &
            "*** Reading Monthly Climatological BC from file " &
            ,trim(specified_bc_file)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        if ( atm_bc_from_netcdf ) then
          !--- open the file containing GT anomaly (if file exists)
          if ( gtano_from_netcdf ) then
            call ncdf_open(fid_gtano, gtano_file_name, mode="read")
          endif

          !--- open the netcdf file containing monthly climatologies
          call ncdf_open(fid_clim, clim_file_name, mode="read")
        endif

        if ( atm_bc_from_cccma ) then
          !--- open the file containing monthly averaged GT,SIC,SICN
          call cc_open(nf21,specified_bc_file,remote,verbose)
        endif

        !--- The current day of the year is iday
        iday=day

        !--- Set delt and gmt such that phsil8 will interpolate to
        !--- the middle of the current day
        delt=0.0
        gmt=43200.0  !--- seconds in the current day

        !--- Define dimension sizes and array limits for PHSIL8
        ilg=lon*lat
        il1=1
        il2=ilg

        !--- Determine mday, the first mid month value after the current day
        !--- and define mon1 and mon2, the two months to interpolate between
        !--- Determine weights for linear interpolation between mid month values
        w1=0.0
        w1=0.0
        mday=0
        mon1=0
        mon2=0
        do k=1,12
          if (k.eq.1) then
            mon1=12
            mon2=1
            if ( (iday.ge.mmd(mon1) .and. iday.le.365) .or. &
                 (iday .ge. 1 .and. iday.lt.mmd(mon2)) ) then
              mday=mmd(mon2)
              dd=(gmt+delt)/86400.0
              den=365 - mmd(mon1) + mmd(mon2)
              if (iday.ge.mmd(mon1) .and. iday.le.365) then
                !--- Last half of Dec
                w2=(float(iday-mmd(mon1))+dd)/den
              else
                !--- First half of Jan
                w2=(float(365-mmd(mon1)+iday)+dd)/den
              endif
              w1=1.0-w2
              exit
            endif
          else
            mon1=k-1
            mon2=k
            if ( iday.ge.mmd(mon1) .and. iday.lt.mmd(mon2)) then
              mday=mmd(mon2)
              dd=(gmt+delt)/86400.0
              den=mmd(mon2)-mmd(mon1)
              w2=(float(iday-mmd(mon1))+dd)/den
              w1=1.0-w2
              exit
            endif
          endif
        enddo
        if (mday.eq.0) call xit('READ_SPEC_BC',-11)

        !--- Determine the current month
        if (iday.ge.som(mon1) .and. iday.le.eom(mon1)) then
          curr_mon=mon1
        else if (iday.ge.som(mon2) .and. iday.le.eom(mon2)) then
          curr_mon=mon2
        else
          curr_mon=0
        endif

        if (verbose.gt.10) then
          strng=' '
          write(strng,'(a,i4,a,i3,a,3i5,a,2g16.8)') &
            ' year=',year,'  day=',day,'  mon1,mon2,mday: ', &
            mon1,mon2,mday,'  w1,w2: ',w1,w2
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name , caller_line, &
                       year=year, mon=curr_mon)
        endif

        !--- Get data from the monthly climatology file for GT, SIC and SICN
        !--- for the two months required for interpolation
        !--- gtrom, sicrom, sicnrom will contain previous mid month data
        !--- gtrol, sicrol, sicnrol will contain next mid month data (ie for mday)

        !--- gtrom --- previous mid month value
        if ( atm_bc_from_cccma ) then
          rewind(nf21)
          call getfld2(nf21,gtrom,-1,mmd(mon1), &
                       nc4to8("  GT"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-12)
          if (verbose.gt.1) then
            gavg = sum(gtrom)/size(gtrom)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(1, mon1)
          call ncdf_read_var(fid_clim, vid_gt, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, gtrom, "gtrom")
          endif
          nx1 = size(gtrom,dim=1)
          ny1 = size(gtrom,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            gtrom = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in gtrom"
            write(6,*)"    shape(gtrom)=",nx1,ny1
            write(6,*)"    shape(data) =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- gtrol --- next mid month value
        if ( atm_bc_from_cccma ) then
          rewind(nf21)
          call getfld2(nf21,gtrol,-1,mmd(mon2), &
                       nc4to8("  GT"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-13)
          if (verbose.gt.1) then
            gavg = sum(gtrol)/size(gtrol)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(1, mon2)
          call ncdf_read_var(fid_clim, vid_gt, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, gtrol, "gtrol")
          endif
          nx1 = size(gtrol,dim=1)
          ny1 = size(gtrol,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            gtrol = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in gtrol"
            write(6,*)"    shape(gtrol)=",nx1,ny1
            write(6,*)"    shape(data) =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- sicrom --- previous mid month value
        if ( atm_bc_from_cccma ) then
          rewind(nf21)
          call getfld2(nf21,sicrom,-1,mmd(mon1), &
                       nc4to8(" SIC"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-14)
          if (verbose.gt.1) then
            gavg = sum(sicrom)/size(sicrom)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(1, mon1)
          call ncdf_read_var(fid_clim, vid_sic, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, sicrom, "sicrom")
          endif
          nx1 = size(sicrom,dim=1)
          ny1 = size(sicrom,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            sicrom = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in sicrom"
            write(6,*)"    shape(sicrom)=",nx1,ny1
            write(6,*)"    shape(data)  =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- sicrol --- next mid month value
        if ( atm_bc_from_cccma ) then
          rewind(nf21)
          call getfld2(nf21,sicrol,-1,mmd(mon2), &
                       nc4to8(" SIC"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-15)
          if (verbose.gt.1) then
            gavg = sum(sicrol)/size(sicrol)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(1, mon2)
          call ncdf_read_var(fid_clim, vid_sic, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, sicrol, "sicrol")
          endif
          nx1 = size(sicrol,dim=1)
          ny1 = size(sicrol,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            sicrol = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in sicrol"
            write(6,*)"    shape(sicrol)=",nx1,ny1
            write(6,*)"    shape(data)  =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- sicnrom --- previous mid month value
        if ( atm_bc_from_cccma ) then
          rewind(nf21)
          call getfld2(nf21,sicnrom,-1,mmd(mon1), &
                       nc4to8("SICN"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-16)
          if (verbose.gt.1) then
            gavg = sum(sicnrom)/size(sicnrom)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(1, mon1)
          call ncdf_read_var(fid_clim, vid_sicn, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, sicnrom, "sicnrom")
          endif
          nx1 = size(sicnrom,dim=1)
          ny1 = size(sicnrom,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            sicnrom = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in sicnrom"
            write(6,*)"    shape(sicnrom)=",nx1,ny1
            write(6,*)"    shape(data)   =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- sicnrol --- next mid month value
        if ( atm_bc_from_cccma ) then
          rewind(nf21)
          call getfld2(nf21,sicnrol,-1,mmd(mon2), &
                       nc4to8("SICN"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-17)
          if (verbose.gt.1) then
            gavg = sum(sicnrol)/size(sicnrol)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(1, mon2)
          call ncdf_read_var(fid_clim, vid_sicn, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, sicnrol, "sicnrol")
          endif
          nx1 = size(sicnrol,dim=1)
          ny1 = size(sicnrol,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            sicnrol = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in sicnrol"
            write(6,*)"    shape(sicnrol)=",nx1,ny1
            write(6,*)"    shape(data)   =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- Interpolate GT to the current day so that input
        !--- points are defined with reasonable values
        !--- Ocean points are overwritten in PHSIL8
        gt(:,:) = w1*gtrom(:,:) + w2*gtrol(:,:)

        !--- Zero out SIC and SICN so that they will be zero over land
        !--- Ocean point are overwritten in PHSIL8
        sic(:,:)  = 0.0
        sicn(:,:) = 0.0

        !--- read GT anomaly to be added to GT after interpolation and clipping in PHSIL
        !--- gtanom --- next mid month value
        if ( gtano_from_netcdf ) then
          time_idx_gtano = min(time_slice_gtano(1850, mon1),ntime_gtano)
          call ncdf_read_var(fid_gtano, vid_gtano, vals_2dr8, nt=time_idx_gtano)
          nx1 = size(gtano,dim=1)
          ny1 = size(gtano,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            gtanom = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in gtano"
            write(6,*)"    shape(gtano)  =",nx1,ny1
            write(6,*)"    shape(data)   =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        else
          gtanom(:,:) = 0.0
        endif

        !--- gtanol --- next mid month value
        if ( gtano_from_netcdf ) then
          time_idx_gtano = min(time_slice_gtano(1850, mon2),ntime_gtano)
          call ncdf_read_var(fid_gtano, vid_gtano, vals_2dr8, nt=time_idx_gtano)
          nx1 = size(gtano,dim=1)
          ny1 = size(gtano,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            gtanol = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in gtano"
            write(6,*)"    shape(gtano)  =",nx1,ny1
            write(6,*)"    shape(data)   =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        else
          gtanol(:,:) = 0.0
        endif

        !--- Interpolate GTano to the current day
        if (ntime_gtano == 1) then
          gtano(:,:) = gtanom(:,:)
        else
          gtano(:,:) = w1*gtanom(:,:) + w2*gtanol(:,:)
        endif

        if (verbose.gt.5) then
          gavg = sum(gt)/size(gt)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' BEFORE PHSIL: year=',year,'  day=',day,'  avg   GT=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
          gavg = sum(sic)/size(sic)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' BEFORE PHSIL: year=',year,'  day=',day,'  avg  SIC=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
          gavg = sum(sicn)/size(sicn)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' BEFORE PHSIL: year=',year,'  day=',day,'  avg SICN=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        call phsil8 (  gt,   gtrol,   gtrom, gtano, &
                      sic,  sicrol,  sicrom, &
                     sicn, sicnrol, sicnrom, &
                     flnd,                   &
                     ilg,il1,il2,delt,gmt,iday,mday,sicn_crt)

        if (verbose.gt.5) then
          gavg = sum(gt)/size(gt)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day,'  avg   GT=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line, &
                       year=year, mon=curr_mon)
          gavg = sum(sic)/size(sic)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day,'  avg  SIC=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line, &
                       year=year, mon=curr_mon)
          gavg = sum(sicn)/size(sicn)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day,'  avg SICN=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line, &
                       year=year, mon=curr_mon)
        endif

        if (verbose.gt.5) then
          strng=' '
          write(strng,'(a,i4,a,i3,a,1p4g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day, &
           '  flnd,gt,sic,sicn at (34,9): ', &
           flnd(34,9),gt(34,9),sic(34,9),sicn(34,9)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line, &
                       year=year, mon=curr_mon)
          strng=' '
          write(strng,'(a,i4,a,i3,a,1p4g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day, &
           '  gtrom,sicrom,sicnrom at (34,9): ', &
           gtrom(34,9),sicrom(34,9),sicnrom(34,9)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line, &
                       year=year, mon=curr_mon)
          strng=' '
          write(strng,'(a,i4,a,i3,a,1p4g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day, &
           '  gtrol,sicrol,sicnrol at (34,9): ', &
           gtrol(34,9),sicrol(34,9),sicnrol(34,9)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line, &
                       year=year, mon=curr_mon)
        endif

        if ( atm_bc_from_netcdf ) then
          call ncdf_close(fid_clim)
          if (gtano_from_netcdf) then
            call ncdf_close(fid_gtano)
          endif
        endif

        if (trim(adjustl(specified_bc_file)) .ne. "AN") then
          !--- Do not close the "AN" file but close all others
          !--- Close the file containing specified boundary conditions
          !--- If this file was "access"ed by cc_open above then the link
          !--- (or copy) will remain in the cwd and it will not need to be
          !--- "access"ed again (as long as the cwd exists).
          if ( atm_bc_from_cccma ) then
            close(nf21)
          endif
        endif

      endif  !--- Monthly Climatologies

      !.......................................................................

      if (specified_bc .eq. 2) then
        !-----------------------------------------------------
        !--- Daily averaged climatological boundary conditions
        !-----------------------------------------------------

        if (verbose.gt.0) then
          write(strng,'(2a)') &
            "*** Reading Daily Climatological BC from file " &
            ,trim(specified_bc_file)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        !--- open the file containing monthly averaged GT,SIC,SICN
        if (len_trim(specified_bc_file) .le. 0) then
          !--- the file name must be supplied by the user
          call xit("READ_SPEC_BC",-21)
        endif
        call cc_open(nf21,specified_bc_file,remote,verbose)

        !--- For climatological data the ibuf2 value is simply
        !--- the day of the year
        timestmp = day

        !--- Read daily average GT, SIC, SICN from a file
        rewind(nf21)
        call getfld2(nf21,gt,-1,timestmp,nc4to8("  GT"), &
                     1,ibuf,maxx,ok)
        if (.not. ok) call xit("READ_SPEC_BC",-22)

        if (verbose.gt.1) then
          strng=' '
          write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8)')ibuf(1:8)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        rewind(nf21)
        call getfld2(nf21,sic,-1,timestmp,nc4to8(" SIC"), &
                     1,ibuf,maxx,ok)
        if (.not. ok) call xit("READ_SPEC_BC",-23)

        if (verbose.gt.1) then
          strng=' '
          write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8)')ibuf(1:8)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        rewind(nf21)
        call getfld2(nf21,sicn,-1,timestmp,nc4to8("SICN"), &
                     1,ibuf,maxx,ok)
        if (.not. ok) call xit("READ_SPEC_BC",-24)

        if (verbose.gt.1) then
          strng=' '
          write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8)')ibuf(1:8)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        !--- Close the file containing specified boundary conditions
        !--- If this file was "access"ed by cc_open above then the link
        !--- (or copy) will remain in the cwd and it will not need to be
        !--- "access"ed again (as long as the cwd exists).
        close(nf21)

        do j=1,lat
          do i=1,lon
            if(flnd(i,j).lt.1. .and. sicn(i,j).eq.0.) then
              !--- Force SNO to be zero over open water
              !--- This will avoid sporadic floating point errors in
              !--- class resulting from zero divides by rhosno
              sno(i,j) = 0.0
              !--- Ensure that ocean surface temperature is not below the
              !--- the freezing point of sea water (GTFSW) at any ocean point
              gt(i,j) = max( gtfsw, gt(i,j) )
            endif
          enddo
        enddo

      endif  !--- Daily Climatologies

      !.......................................................................

      if (specified_bc .eq. 3) then
        !-------------------------------------------------------
        !--- Monthly multi-year boundary conditions
        !-------------------------------------------------------

        if (verbose.gt.0) then
          write(strng,'(2a)') &
            "*** Reading monthly multi-year SST/SIC/SICN from file " &
            ,trim(specified_bc_file)
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        if (len_trim(specified_bc_file) .le. 0) then
          !--- the file name must be supplied by the user
          call xit("READ_SPEC_BC",-31)
        endif

        if ( atm_bc_from_netcdf ) then
          !--- open the file containing GT anomaly (if file exists)
          if ( gtano_from_netcdf ) then
            call ncdf_open(fid_gtano, gtano_file_name, mode="read")
          endif

          !--- open the file containing monthly averaged multi-year SST
          call ncdf_open(fid_gt, gt_file_name, mode="read")

          !--- open the file containing monthly averaged multi-year SIC
          call ncdf_open(fid_sic, sic_file_name, mode="read")

          !--- open the file containing monthly averaged multi-year SICN
          call ncdf_open(fid_sicn, sicn_file_name, mode="read")
        endif

        if ( atm_bc_from_cccma ) then
          !--- open the file containing monthly averaged multi-year SST
          strng=trim(adjustl(specified_bc_file))//"_gt"
          call cc_open(nf21,strng,remote,verbose)
          rewind nf21

          !--- open the file containing monthly averaged multi-year SIC
          strng=trim(adjustl(specified_bc_file))//"_sic"
          call cc_open(nf22,strng,remote,verbose)
          rewind nf22

          !--- open the file containing monthly averaged multi-year SICN
          strng=trim(adjustl(specified_bc_file))//"_sicn"
          call cc_open(nf23,strng,remote,verbose)
          rewind nf23
        endif

        !--- The current day of the year is iday
        iday=day

        !--- The current year is iyear
        !--- specified_bc_year_offset is defined in com_cpl
        if ( specified_bc_years_in_loop > 0 ) then
!TODO--- Define year_offset such that the multi year boundary conditions will repeat
!TODO--- every specified_bc_years_in_loop years
write(6,*)"read_spec_bc: specified_bc_years_in_loop > 0 is not yet supported."
call xit("READ_SPEC_BC",-90)
          year_offset = (year - specified_bc_start_year) + specified_bc_year_offset
        else
          year_offset = specified_bc_year_offset
        endif
        iyear = year + year_offset - 1

        !--- Set delt and gmt such that phsil8 will interpolate to
        !--- the middle of the current day
        delt=0.0
        gmt=43200.0  !--- seconds in the current day

        !--- Define dimension sizes and array limits for PHSIL8
        ilg=lon*lat
        il1=1
        il2=ilg

write(6,*)"read_spec_bc: iyear,iday=",iyear,iday
call flush(6)

        !--- Determine mday, the first mid month value after the current day
        !--- and define mon1 and mon2, the two months to interpolate between
        !--- Determine weights for linear interpolation between mid month values
        w1=0.0
        w1=0.0
        mday=0
        mon1=0
        mon2=0
        do k=1,12
          if (k.eq.1) then
            mon1=12
            mon2=1
            if ( (iday.ge.mmd(mon1) .and. iday.le.365) .or. &
                 (iday .ge. 1 .and. iday.lt.mmd(mon2)) ) then
              mday=mmd(mon2)
              dd=(gmt+delt)/86400.0
              den=365 - mmd(mon1) + mmd(mon2)
              if (iday.ge.mmd(mon1) .and. iday.le.365) then
                !--- Last half of Dec
                w2=(float(iday-mmd(mon1))+dd)/den
              else
                !--- First half of Jan
                w2=(float(365-mmd(mon1)+iday)+dd)/den
              endif
              w1=1.0-w2
              exit
            endif
          else
            mon1=k-1
            mon2=k
            if ( iday.ge.mmd(mon1) .and. iday.lt.mmd(mon2)) then
              mday=mmd(mon2)
              dd=(gmt+delt)/86400.0
              den=mmd(mon2)-mmd(mon1)
              w2=(float(iday-mmd(mon1))+dd)/den
              w1=1.0-w2
              exit
            endif
          endif
        enddo
        if (mday.eq.0) call xit('READ_SPEC_BC',-32)

        !--- Determine the ibuf(2) values of the *rom and *rol records to be
        !--- read from the muti year GT, SIC and SICN data files
        if (mon1.eq.12.and.day.le.15) then
          !--- The first month is Dec of the previous year
          ib2rom = (iyear-1)*1000000 + mon1*10000 + 100
          year1 = iyear-1
        else
          !--- The first month is mon1 of the current year
          ib2rom = iyear*1000000 + mon1*10000 + 100
          year1 = iyear
        endif
        if (mon1.eq.12.and.day.ge.350) then
          !--- The second month is mon2 of the next year
          ib2rol = (iyear+1)*1000000 + mon2*10000 + 100
          year2 = iyear+1
        else
          !--- The second month is mon2 of the current year
          ib2rol = iyear*1000000 + mon2*10000 + 100
          year2 = iyear
        endif

        if (verbose.gt.10) then
          strng=' '
          write(strng,'(a,i4,a,i3,a,3i5,a,2g16.8)') &
            ' year=',year,'  day=',day,'  mon1,mon2,mday: ', &
            mon1,mon2,mday,'  w1,w2: ',w1,w2
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        !--- Get data from the monthly climatology file for GT, SIC and SICN
        !--- for the two months required for interpolation
        !--- gtrom, sicrom, sicnrom will contain previous mid month data
        !--- gtrol, sicrol, sicnrol will contain next mid month data (ie for mday)

        !--- gtrom --- previous mid month data
        if ( atm_bc_from_cccma ) then
          rewind(nf21)
          call getfld2(nf21,gtrom,-1,ib2rom, &
                       nc4to8("  GT"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-33)
          if (verbose.gt.1) then
            gavg = sum(gtrom)/size(gtrom)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(year1, mon1)
          call ncdf_read_var(fid_gt, vid_gt, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, gtrom, "gtrom")
          endif
          nx1 = size(gtrom,dim=1)
          ny1 = size(gtrom,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            gtrom = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in gtrom"
            write(6,*)"    shape(gtrom)=",nx1,ny1
            write(6,*)"    shape(data) =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- gtrol --- next mid month data
        if ( atm_bc_from_cccma ) then
          rewind(nf21)
          call getfld2(nf21,gtrol,-1,ib2rol, &
                       nc4to8("  GT"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-34)
          if (verbose.gt.1) then
            gavg = sum(gtrol)/size(gtrol)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(year2, mon2)
          call ncdf_read_var(fid_gt, vid_gt, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, gtrol, "gtrol")
          endif
          nx1 = size(gtrol,dim=1)
          ny1 = size(gtrol,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            gtrol = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in gtrol"
            write(6,*)"    shape(gtrol)=",nx1,ny1
            write(6,*)"    shape(data) =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- sicrom --- previous mid month data
        if ( atm_bc_from_cccma ) then
          rewind(nf22)
          call getfld2(nf22,sicrom,-1,ib2rom, &
                       nc4to8(" SIC"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-35)
          if (verbose.gt.1) then
            gavg = sum(sicrom)/size(sicrom)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(year1, mon1)
          call ncdf_read_var(fid_sic, vid_sic, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, sicrom, "sicrom")
          endif
          nx1 = size(sicrom,dim=1)
          ny1 = size(sicrom,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            sicrom = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in sicrom"
            write(6,*)"    shape(sicrom)=",nx1,ny1
            write(6,*)"    shape(data)  =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- sicrol --- next mid month data
        if ( atm_bc_from_cccma ) then
          rewind(nf22)
          call getfld2(nf22,sicrol,-1,ib2rol, &
                       nc4to8(" SIC"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-36)
          if (verbose.gt.1) then
            gavg = sum(sicrol)/size(sicrol)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(year2, mon2)
          call ncdf_read_var(fid_sic, vid_sic, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, sicrol, "sicrol")
          endif
          nx1 = size(sicrol,dim=1)
          ny1 = size(sicrol,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            sicrol = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in sicrol"
            write(6,*)"    shape(sicrol)=",nx1,ny1
            write(6,*)"    shape(data)  =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- sicnrom --- previous mid month data
        if ( atm_bc_from_cccma ) then
          rewind(nf23)
          call getfld2(nf23,sicnrom,-1,ib2rom, &
                       nc4to8("SICN"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-37)
          if (verbose.gt.1) then
            gavg = sum(sicnrom)/size(sicnrom)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(year1, mon1)
          call ncdf_read_var(fid_sicn, vid_sicn, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_netcdf .and. atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, sicnrom, "sicnrom")
          endif
          nx1 = size(sicnrom,dim=1)
          ny1 = size(sicnrom,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            sicnrom = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in sicnrom"
            write(6,*)"    shape(sicnrom)=",nx1,ny1
            write(6,*)"    shape(data)   =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- sicnrol --- next mid month data
        if ( atm_bc_from_cccma ) then
          rewind(nf23)
          call getfld2(nf23,sicnrol,-1,ib2rol, &
                       nc4to8("SICN"),1,ibuf,maxx,ok)
          if (.not. ok) call xit('READ_SPEC_BC',-38)
          if (verbose.gt.1) then
            gavg = sum(sicnrol)/size(sicnrol)
            strng=' '
            write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g16.8)') &
              ibuf(1:8),'  avg=',gavg
            caller_line      = __LINE__
            call dbg_out(strng, caller_file_name, caller_line)
          endif
        endif

        if ( atm_bc_from_netcdf ) then
          time_idx = time_slice(year2, mon2)
          call ncdf_read_var(fid_sicn, vid_sicn, vals_2dr8, nt=time_idx)
          if ( atm_bc_from_cccma ) then
            call compare_2d_r8(vals_2dr8, sicnrol, "sicnrol")
          endif
          nx1 = size(sicnrol,dim=1)
          ny1 = size(sicnrol,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            sicnrol = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in sicnrol"
            write(6,*)"    shape(sicnrol)=",nx1,ny1
            write(6,*)"    shape(data)   =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        endif

        !--- Interpolate GT to the current day so that input
        !--- points are defined with reasonable values
        !--- Ocean points are overwritten in PHSIL8
        gt(:,:) = w1*gtrom(:,:) + w2*gtrol(:,:)

        !--- Zero out SIC and SICN so that they will be zero over land
        !--- Ocean point are overwritten in PHSIL8
        sic(:,:)  = 0.0
        sicn(:,:) = 0.0

        !--- read GT anomaly to be added to GT after interpolation and clipping in PHSIL
        !--- gtanom --- previous mid month value
        if ( gtano_from_netcdf ) then
          time_idx_gtano = min(time_slice_gtano(1850, mon1),ntime_gtano)
          call ncdf_read_var(fid_gtano, vid_gtano, vals_2dr8, nt=time_idx_gtano)
          nx1 = size(gtano,dim=1)
          ny1 = size(gtano,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            gtanom = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in gtano"
            write(6,*)"    shape(gtano)  =",nx1,ny1
            write(6,*)"    shape(data)   =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        else
          gtanom(:,:) = 0.0
        endif

        !--- gtanol --- next mid month value
        if ( gtano_from_netcdf ) then
          time_idx_gtano = min(time_slice_gtano(1850, mon2),ntime_gtano)
          call ncdf_read_var(fid_gtano, vid_gtano, vals_2dr8, nt=time_idx_gtano)
          nx1 = size(gtano,dim=1)
          ny1 = size(gtano,dim=2)
          nx2 = size(vals_2dr8,dim=1)
          ny2 = size(vals_2dr8,dim=2)
          if ( nx1 == nx2 .and. ny1 == ny2 ) then
            gtanol = vals_2dr8
          else
            write(6,*)"read_spec_bc: shape mismatch in gtano"
            write(6,*)"    shape(gtano)  =",nx1,ny1
            write(6,*)"    shape(data)   =",nx2,ny2
            call xit("read_spec_bc",-91)
          endif
          deallocate(vals_2dr8)
        else
          gtanol(:,:) = 0.0
        endif

        !--- Interpolate GTano to the current day
        if (ntime_gtano == 1) then
          gtano(:,:) = gtanom(:,:)
        else
          gtano(:,:) = w1*gtanom(:,:) + w2*gtanol(:,:)
        endif
        
        if (verbose.gt.5) then
          gavg = sum(gt)/size(gt)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' BEFORE PHSIL: year=',year,'  day=',day,'  avg   GT=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
          gavg = sum(sic)/size(sic)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' BEFORE PHSIL: year=',year,'  day=',day,'  avg  SIC=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
          gavg = sum(sicn)/size(sicn)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' BEFORE PHSIL: year=',year,'  day=',day,'  avg SICN=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        call phsil8 (  gt,   gtrol,   gtrom, gtano, &
                      sic,  sicrol,  sicrom, &
                     sicn, sicnrol, sicnrom, &
                     flnd,                   &
                     ilg,il1,il2,delt,gmt,iday,mday,sicn_crt)

        if (verbose.gt.5) then
          gavg = sum(gt)/size(gt)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day,'  avg   GT=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
          gavg = sum(sic)/size(sic)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day,'  avg  SIC=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
          gavg = sum(sicn)/size(sicn)
          strng=' '
          write(strng,'(a,i4,a,i3,a,g16.8)') &
           ' AFTER PHSIL: year=',year,'  day=',day,'  avg SICN=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        !--- Close the files containing specified boundary conditions
        !--- If these files were "access"ed by cc_open above then the link
        !--- (or copy) will remain in the cwd and they will not need to be
        !--- "access"ed again (as long as the cwd exists).
        if ( atm_bc_from_netcdf ) then
          call ncdf_close(fid_gt)
          call ncdf_close(fid_sic)
          call ncdf_close(fid_sicn)
          if (gtano_from_netcdf) then
            call ncdf_close(fid_gtano)
          endif
        endif

        if ( atm_bc_from_cccma ) then
          close(nf21)
          close(nf22)
          close(nf23)
        endif

      endif  !--- Monthly multi-year boundary conditions

      !.......................................................................

      if (specified_bc .eq. 4) then
        !-------------------------------------------------------
        !--- Daily multi-year boundary conditions
        !-------------------------------------------------------

        !--- Given the day of the year determine the corresponding
        !--- month and day of the month
        month=0
        dom=0
        do i=1,12
          if (som(i).le.day .and. eom(i).ge.day) then
            month = i
            dom   = 1 + day - som(i)
            exit
          endif
        enddo
        if (month.lt.1 .or. month.gt.12) call xit("READ_SPEC_BC",-41)
        if (dom.lt.1 .or. dom.gt.31)     call xit("READ_SPEC_BC",-42)

        !--- Determine the timestamp (ibuf2 value) of the records to
        !--- be read from the daily data files
        timestmp = year*1000000 + month*10000 + dom*100

        if (verbose.gt.0) then
          write(strng,'(3a,i12)') &
            "*** Reading daily multi-year SST/SIC/SICN from file " &
            ,trim(specified_bc_file) &
            ,'   ibuf2=',timestmp
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
          write(strng,'(a,3g12.4)')" day,month,dom: ",day,month,dom
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        if (len_trim(specified_bc_file) .le. 0) then
          !--- the file name must be supplied by the user
          call xit("READ_SPEC_BC",-43)
        endif

        !--- open the file containing daily averaged multi-year GT
        strng=trim(adjustl(specified_bc_file))//"_gt"
        call cc_open(nf21,strng,remote,verbose)
        rewind nf21

        call getfld2(nf21,gt,-1,timestmp, &
                     nc4to8("  GT"),1,ibuf,maxx,ok)
        if (.not. ok) call xit("READ_SPEC_BC",-44)
        close(nf21)

        if (verbose.gt.1) then
          gavg = sum(GT)/size(GT)
          strng=' '
          write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g12.4)') &
            ibuf(1:8),'  avg=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        !--- open the file containing daily averaged multi-year SIC
        strng=trim(adjustl(specified_bc_file))//"_sic"
        call cc_open(nf21,strng,remote,verbose)
        rewind nf21

        call getfld2(nf21,sic,-1,timestmp, &
                     nc4to8(" SIC"),1,ibuf,maxx,ok)
        if (.not. ok) call xit("READ_SPEC_BC",-45)
        close(nf21)

        if (verbose.gt.1) then
          gavg = sum(SIC)/size(SIC)
          strng=' '
          write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g12.4)') &
            ibuf(1:8),'  avg=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        !--- open the file containing daily averaged multi-year SICN
        strng=trim(adjustl(specified_bc_file))//"_sicn"
        call cc_open(nf21,strng,remote,verbose)
        rewind nf21

        call getfld2(nf21,sicn,-1,timestmp, &
                     nc4to8("SICN"),1,ibuf,maxx,ok)
        if (.not. ok) call xit("READ_SPEC_BC",-46)
        close(nf21)

        if (verbose.gt.1) then
          gavg = sum(SICN)/size(SICN)
          strng=' '
          write(strng,'(2x,a4,1x,i10,1x,a4,1x,5i8,a,g12.4)') &
            ibuf(1:8),'  avg=',gavg
          caller_line      = __LINE__
          call dbg_out(strng, caller_file_name, caller_line)
        endif

        do j=1,lat
          do i=1,lon
            if(flnd(i,j).lt.1. .and. sicn(i,j).eq.0.) then
              !--- Force SNO to be zero over open water
              !--- This will avoid sporadic floating point errors in
              !--- class resulting from zero divides by rhosno
              sno(i,j) = 0.0
              !--- Ensure that ocean surface temperature is not below the
              !--- the freezing point of sea water (GTFSW) at any ocean point
              gt(i,j) = max( gtfsw, gt(i,j) )
            endif
          enddo
        enddo

      endif  !--- Daily multi-year boundary conditions

      !.......................................................................
      !--- Force SNO to be zero over open water
      do j=1,lat
        do i=1,lon
          if ( (flnd(i,j) .lt. 1.) .and. (sicn(i,j) .eq. 0.) ) then
            !--- No ice, open water ocean grid cell
            sno(i,j) = 0.0
          endif
        enddo
      enddo

      check_ice = .false.
      if (check_ice) then
        !--- Adjust values of SIC so they are consistent with those of SICN
        !--- Limit ice thickness to 6 meters
        do i=1,lon
          do j=1,lat

            if (flnd(i,j) .lt. 1.) then
              !--- Ensure both SIC and SICN are zero over land
              sic (i,j)=0.0
              sicn(i,j)=0.0 !--- also make inflated values consistent
            endif

            if (sicn(i,j) .le. 0.0) then
              !--- Ensure SIC=0 when SICN=0
              sic(i,j)=0.0
            endif

            if (sicn(i,j) .lt. sicn_crt) then
              !--- When SICN < SICN_CRT set SIC = a number between 0 and 30
              sic(i,j) = min(sic(i,j),30.0)
            endif

            if (sicn(i,j) .eq. sicn_crt) then
              !--- When SICN = SICN_CRT set SIC = SIC_CRT = 45
              sic(i,j) = 45.0
            endif

            if (sicn(i,j) .gt. sicn_crt) then
              !--- When SICN > SICN_CRT set SIC = a number greater than SIC_CRT
              sic(i,j) = max(sic(i,j),45.0)

              !--- Ensure SIC will produce a maximum ice thickness of 6m
              !--- ice thickness:       h = SIC/(rho_ice*SICN)
              !--- ice density:   rho_ice = 913 kg/m**3
              sicmax=6.0*913.0*sicn(i,j)
              sic(i,j) = min(sic(i,j),sicmax)
            endif

          enddo
        enddo
      endif

    end subroutine read_spec_bc

    !***********************************************************************
    ! Write coupler debugging output to a text file
    !***********************************************************************
    subroutine dbg_out(msg, caller, line, year, mon)
      character*(*) :: msg, caller
      integer :: line
      integer, optional :: year, mon

      !--- Local
      character(256),save :: runid=' ', fout=' '
      integer, save :: NFDBG = 557
      logical, save :: first_call=.true.

      character(256) :: usr_home, fname
      character(10)  :: curr_date=' ', curr_time=' '

      !--- Set runid
#ifdef _RUNID_
      !--- The user has supplied a value for _RUNID_
#ifdef CrayFTN
#define QUOTEMAC(x) "x"
#else
#define QUOTEMAC(x) #x
#endif
#define QUOTEMACVAL(x) QUOTEMAC(x)
      runid= QUOTEMACVAL(_RUNID_)
#else
      runid=' '
#endif

      if (first_call) then
        first_call = .false.

        !--- Determine the output file name
        if ( present(year) .and. present(mon) ) then
          write(fname,'(a,i4.4,a,i2.2)')'dbg_out_cpl_',year,'_m',mon
        else if ( present(year) ) then
          write(fname,'(a,i4.4)')'dbg_out_cpl_',year
        else if (present(mon) ) then
          write(fname,'(a,i2.2)')'dbg_out_cpl_m',mon
        else
          fname='dbg_out_cpl'
        endif
        if (len_trim(runid).gt.0) then
          fname=trim(fname)//'_'//trim(adjustl(runid))
        endif

!!$        usr_home=' '
!!$        call getenv('HOME', usr_home)
!!$        if (len_trim(usr_home).gt.0) then
!!$          !--- If we can get the value of the environment variable HOME
!!$          !--- then create the debug output file relative to the
!!$          !--- invoking users home dir
!!$          fout=trim(adjustl(usr_home))//'/tmp/'//trim(adjustl(fname))
!!$        else
          !--- Otherwise the debug output file goes into the cwd
          fout=trim(adjustl(fname))
!!$        endif

        !--- Create the output file, overwriting any existing file
        open(NFDBG,file=trim(fout),form='formatted',position='REWIND')
        call date_and_time(curr_date, curr_time)
        write(NFDBG,'(4a)')" DATE: ",curr_date,"   TIME: ",curr_time
        close(NFDBG)

      endif

      open(NFDBG,file=trim(fout),form='formatted',position='APPEND')
      if (len_trim(msg).gt.0) then
        write(NFDBG,'(2a,i6,1x,a)')trim(caller),':',line,trim(msg)
      else
        write(NFDBG,'(2a,i6)')trim(caller),':',line
      endif
      close(NFDBG)

    end subroutine dbg_out

    !***********************************************************************
    !--- Extract info about time units and time origin from the input string
    !***********************************************************************
    subroutine parse_time_units(time_units_in)

      implicit none

      character(*) :: time_units_in

      !--- local
      character(len=64) :: units, strng, word(10), numb(10)
      integer :: idx1, idx2, idx3, verbose=1
      logical :: ok
      integer :: year, mon, day, hour, min, sec, nword, nnumb

      !--- The input time_units_in must be a string of the form
      !---    U since Y-M-D H:M:S
      !--- where U is one of "days", "hours", "minutes" or "seconds"
      !--- certain abbreviations of these units are allowed

      if (len_trim(time_units_in).le.0) then
        write(6,*)"parse_time_units: Empty input string"
        call xit("PARSE_TIME_UNITS",-1)
      endif

      !--- Initialize variables that may be set below
      units=" "
      year=0
      mon=1
      day=1
      hour=0
      min=0
      sec=0

      !--- Make everything lower case to simplify the comparisons done below
      strng = " "
      call tolower(time_units_in,strng)

      !--- Remove any leading white space in time_units_in
      strng = adjustl(strng)

      !--- Separate the time_units_in string on white space into an array of words
      !--- call str2list(time_units_in,64," ",word,nword,4)
      call split_on_delim(word,nword,strng," ")

      if (nword.lt.3) then
        !--- The "U since Y-M-D" portion of the attribute string is mandatory
        !--- The "H:M:S" field is optional
        write(6,'(3a)')"parse_time_units: Invalid time_units string --> ",trim(time_units_in)," <--"
        call xit("PARSE_TIME_UNITS",-2)
      endif

      !--- Process each word to extract time units and origin
      do idx1=1,nword
        select case (idx1)
          case (1)
            !--- This should be the units string
            time_units = " "
            time_units = trim(word(idx1))
            !--- Set time_units explicitly in the following case statement so
            !--- that regardless of the abbreviation, time_units will always
            !--- contain a know string which may be compared later on
            select case (trim(adjustl(time_units)))
              case ("y", "year", "years")
                time_units = "years"
              case ("d", "day", "days")
                time_units = "days"
              case ("h", "hr", "hrs", "hour", "hours")
                time_units = "hours"
              case ("m", "min", "mins", "minute", "minutes")
                time_units = "minutes"
              case ("s", "sec", "secs", "second", "seconds")
                time_units = "seconds"
              case default
                !--- Invalid time units
                write(6,'(3a)')"parse_time_units: Invalid time units --> ", &
                  trim(time_units_in)," <--"
                write(6,'(a)')" The first word must be one of days, hours, minutes, seconds."
                call xit("PARSE_TIME_UNITS",-3)
            end select
          case (2)
            !--- This should be the string "since"
            if (trim(adjustl(word(idx1))) .ne. "since" ) then
              !--- Abort if this is not "since"
              write(6,'(3a)')"parse_time_units: Invalid time units --> ", &
                trim(time_units_in)," <--"
              write(6,'(a)')' The second word must be "since".'
              call xit("PARSE_TIME_UNITS",-4)
            endif
          case (3)
            !--- This should be the date string in the form Y-M-D
            !--- call str2list(word(idx1),64,"-",numb,nnumb,3)
            call split_on_delim(numb,nnumb,word(idx1),"-")
            if (nnumb.ne.3) then
              write(6,'(3a)')"parse_time_units: Invalid time units --> ", &
                trim(time_units_in)," <--"
              write(6,'(3a)')" Invalid date field --> ", &
                trim(word(idx1))," <-- This field must be of the form Y-M-D."
              call xit("PARSE_TIME_UNITS",-5)
            endif
            do idx2=1,nnumb
              strng=" "
              strng=trim(adjustl(numb(idx2)))
              select case (idx2)
                case (1)
                  read(strng,*)year
                  if (year.lt.0) then
                    write(6,*)"parse_time_units: year is out of range. year=",year
                    call xit("PARSE_TIME_UNITS",-6)
                  endif
                case (2)
                  read(strng,*)mon
                  if (mon.lt.1 .or. mon.gt.12) then
                    write(6,*)"parse_time_units: month is out of range. mon=",mon
                    call xit("PARSE_TIME_UNITS",-7)
                  endif
                case (3)
                  read(strng,*)day
                  if (day.lt.1 .or. day.gt.31) then
                    write(6,*)"parse_time_units: day is out of range. day=",day
                    call xit("PARSE_TIME_UNITS",-8)
                  endif
                case default
                  !--- Too many fields
                  write(6,*)"parse_time_units: Too many fields found in date string."
                  write(6,*)"   ",trim(word(idx1))
                  call xit("PARSE_TIME_UNITS",-9)
              end select
            enddo
          case (4)
            !--- This should be the time of day string in the form H:M:S
            !--- call str2list(word(idx1),64,":",numb,nnumb,3)
            call split_on_delim(numb,nnumb,word(idx1),":")
            if (nnumb.ne.3) then
              write(6,'(3a)')"parse_time_units: Invalid time units --> ", &
                trim(time_units_in)," <--"
              write(6,'(3a)')" Invalid time of day field --> ", &
                trim(word(idx1))," <-- This field must be of the form H:M:S."
              call xit("PARSE_TIME_UNITS",-10)
            endif
            do idx2=1,nnumb
              strng=" "
              strng=trim(adjustl(numb(idx2)))
              select case (idx2)
                case (1)
                  read(strng,*)hour
                  if (hour.lt.0 .or. hour.gt.24) then
                    write(6,*)"parse_time_units: hour is out of range. hour=",hour
                    call xit("PARSE_TIME_UNITS",-11)
                  endif
                case (2)
                  read(strng,*)min
                  if (min.lt.0 .or. min.gt.60) then
                    write(6,*)"parse_time_units: minute is out of range. min=",min
                    call xit("PARSE_TIME_UNITS",-12)
                  endif
                case (3)
                  read(strng,*)sec
                  if (sec.lt.0 .or. sec.gt.60) then
                    write(6,*)"parse_time_units: second is out of range. sec=",sec
                    call xit("PARSE_TIME_UNITS",-13)
                  endif
                case default
                  !--- Too many fields
                  write(6,*)"parse_time_units: Too many fields found in time of day string."
                  write(6,*)"   ",trim(word(idx1))
                  call xit("PARSE_TIME_UNITS",-14)
              end select
            enddo
          case default
            !--- Too many words
            write(6,*)"parse_time_units: Too many words found in time_units string."
            write(6,*)"   ",trim(time_units_in)
            call xit("PARSE_TIME_UNITS",-15)
        end select
      enddo

      !--- By this point the parameters year, mon, day, hour, min, sec should have been
      !--- set appropriately. Now assign them to the global reference date/time array
      time_origin_ymd_hms(1) = year
      time_origin_ymd_hms(2) = mon
      time_origin_ymd_hms(3) = day
      time_origin_ymd_hms(4) = hour
      time_origin_ymd_hms(5) = min
      time_origin_ymd_hms(6) = sec

      if ( verbose > 0 ) then
        write(6,*)"parse_time_units: Y M D H M S ",time_origin_ymd_hms
      endif

    end subroutine parse_time_units

    !***********************************************************************
    !--- Find min, max and relative differences between the input arrays
    !***********************************************************************
    subroutine compare_2d_r8(var1, var2, msg)
      real(kind=8), intent(in) :: var1(:,:), var2(:,:)
      character(*), intent(in) :: msg

      !--- Local
      real(kind=8), pointer :: wrk(:,:), wrk2(:,:)
      real(kind=8) :: gsize
      integer :: nx, ny

      nullify(wrk,wrk2)

      if ( size(var1,dim=1) == size(var2,dim=1) .and. &
           size(var1,dim=2) == size(var2,dim=2) ) then
        nx = size(var1,dim=1)
        ny = size(var1,dim=2)
        allocate( wrk(nx,ny), wrk2(nx,ny) )
        wrk = var1 - var2
        write(6,*)trim(msg),": min and max diff ",minval(wrk),maxval(wrk)
        where (wrk < 0.0_8) wrk = -wrk
        wrk2 = var1
        where (wrk2 < 0.0_8) wrk2 = -wrk2
        gsize = sum(wrk2)
        if ( gsize == 0.0_8 ) then
          gsize = nx*ny
          write(6,*)trim(msg),": avg sum of abs diff ",sum(wrk)/gsize
        else
          write(6,*)trim(msg),": relative abs diff ",sum(wrk)/gsize
        endif
        deallocate( wrk, wrk2 )
      else
        write(6,*)trim(msg),": shape mismatch"
        write(6,*)trim(msg),": shape(var1)=",shape(var1)
        write(6,*)trim(msg),": shape(var2)=",shape(var2)
      endif
    end subroutine compare_2d_r8

    !***********************************************************************
    !--- Define PHSIL8 here
    !--- This is the only place that phsil8 should be used
    !***********************************************************************
    SUBROUTINE PHSIL8 (  GTROW,  GTROL,  GTROM, GTANO, &
                        SICROW, SICROL, SICROM, &
                       SICNROW,SICNROL,SICNROM, &
                       FLND,                    &
                       ILG,IL1,IL2,DELT,GMT,IDAY,MDAY,SICN_CRT)
!
!     * FEB 02/19 - S.KHARIN. - ADD GTANO.
!     * AUG 12/18 - M.LAZARE. - REMOVE UNUSED DSICROW.
!     * APR 10/17 - S.KHARIN NEW VERSION FOR GCM18:
!                             - REMOVE GC,MASK AND USE FLND INSTEAD
!                             - SIMPLIFY LOGIC FOR SICN/SIC/GT
!                             - INTRODUCE ADDITIONAL CONSTRAINT ON SIC
!                               FOR SICN<0.15.
!     * APR 23/10 - M.LAZARE. NEW VERSION FOR GCM15I:
!     *                       - SET GTROW=GTROL WHEN END INTERPOLATION
!     *                         TARGET DATE IS REACHED, FOR TRANSITION
!     *                         FROM SEAICE TO WATER, TO AVOID ROUNDOFF
!     *                         ERROR.
!     * APR 02/08 - M.LAZARE. REVISED (COSMETIC) VERSION PHSIL5
!     *                       FOR GCM15F/G/H:
!     *                       - INITIALIZE DSICROW TO ZERO!
!     * JUN 15/06 - M.LAZARE. PREVIOUS VERSION PHSIL5 FOR GCM15F:
!     *                       - "ZERO","ONE" DATA CONSTANTS ADDED AND
!     *                         USED IN CALLS TO INTRINSICS.
!     * NOV 01/04 - M.LAZARE. PEVIOUS VERSION PHSIL4 UP TO AND INCL GCM15E.
!     *                       -  COSMETIC CHANGE: FORCE SICN=0. OVER LAND.
!     * JUN 10/03 - MLAZARE.  CHANGE 1,LON TO IL1,IL2 (PASSED IN).
!     * MAR 17/99 - M.LAZARE. PREVIOUS VERSION PHSIL3.
!
!     * GTROW  = GROUND TEMPERATURE    FOR PREVIOUS TIMESTEP.
!     * GTROL  = GROUND TEMPERATURE    AT NEXT PHYSICS TIME.
!     * GTROM  = GROUND TEMPERATURE    AT PREVIOUS PHYSICS TIME.
!     * SICROW = SEA-ICE MASS          FOR PREVIOUS TIMESTEP.
!     * SICROL = SEA-ICE MASS          AT NEXT PHYSICS TIME.
!     * SICROM = SEA-ICE MASS          AT PREVIOUS PHYSICS TIME.
!     * SICNROW= SEA-ICE CONCENTRATION FOR PREVIOUS TIMESTEP.
!     * SICNROL= SEA-ICE CONCENTRATION AT NEXT PHYSICS TIME.
!     * SICNROM= SEA-ICE CONCENTRATION AT PREVIOUS PHYSICS TIME.
!     * FLND   = FRACTIONAL LAND MASK.
!
!     * ILG    = SIZE OF LONGITUDE SLICE.
!     * IL1    = START OF LONGITUDE SLICE.
!     * IL2    = END OF LONGITUDE SLICE.
!     * DELT   = MODEL TIMESTEP IN SECONDS.
!     * GMT    = NUMBER OF SECONDS IN CURRENT DAY.
!     * IDAY   = CURRENT JULIAN DAY.
!     * DAY    = CURRENT TIME (IDAY PLUS TIME OF DAY).
!     * MDAY   = DATE OF NEXT      MID-MONTH PHYSICS.
!     * MDAYM  = DATE OF PREVIOUS  MID-MONTH PHYSICS.
!     * SICN_CRT=THRESHOLD OF SEA-ICE FRACTION ABOVE WHICH POINT IS
!     *          CONSIDERED TO BE ICE-COVERED.
!
      IMPLICIT REAL (A-H,O-Z), INTEGER (I-N)

      REAL     GTROW(ILG),  GTROL(ILG),  GTROM(ILG), GTANO(ILG)
      REAL    SICROW(ILG), SICROL(ILG), SICROM(ILG)
      REAL   SICNROW(ILG),SICNROL(ILG),SICNROM(ILG)
      REAL   FLND   (ILG)
!
      INTEGER MMD(12)
!
      COMMON /PARAM3/ CSNO,   CPACK, GTFSW,  RKHI, SBC,  SNOMAX
      COMMON /PARAM5/ CONI,   DENI,  SIC_MIN,DFSIC,CONF
!
      DATA  MMD/ 16, 46, 75,106,136,167,197,228,259,289,320,350/
      DATA ZERO,ONE / 0., 1. /
!--------------------------------------------------------------------
      M=0
      DO N=1,12
       IF(MDAY.EQ.MMD(N))M=N
      ENDDO
      IF(M.EQ.0)CALL XIT('PHSIL8',-1)
      MM1=M-1
      IF(MM1.EQ.0)MM1=12
      MDAYM=MMD(MM1)
!
!     * COMPUTE THE TIME DIFFERENCES.
!
      DAY=REAL(IDAY)+(GMT+DELT)/86400.0
      DT1=DAY-REAL(MDAYM)
      DT2=REAL(MDAY)-DAY
      DT3=REAL(MDAY-MDAYM)
      IF(DT1.LT.0.0)DT1=DT1+365.0
      IF(DT2.LT.0.0)DT2=DT2+365.0
      IF(DT3.LT.0.0)DT3=DT3+365.0
!
      DO I=IL1,IL2
!
!     * CASE 1: OVER LAND (NOTE: FORCE SEA-ICE CONCENTRATIONS TO ZERO)
!
       IF(FLND(I).EQ.1.0)THEN
        SICNROW(I)=0.
        SICROW(I)=0.
       ELSE
!
!       * INTERPOLATED SST
!
        GTROW(I)=GTROM(I)+(GTROL(I)-GTROM(I))*DT1/DT3
!
!       * CLIP AT FREEZING POINT AND ADD ANOMALY
!
        GTROW(I)=MAX(GTROW(I),GTFSW)+GTANO(I)
!
!       * INTERPOLATED SICN
!
        SICNROW(I)=SICNROM(I)+(SICNROL(I)-SICNROM(I))*DT1/DT3
!
!       * CLIP AT (0...1)
!
        SICNROW(I)=MAX(ZERO,MIN(SICNROW(I),ONE))
!
!       * INTERPOLATED SIC
!
        SICOLD=SICROW(I)
        SICROW(I)=SICROM(I)+(SICROL(I)-SICROM(I))*DT1/DT3
        IF(SICNROW(I).GT.0.15)THEN
!
!       * RESTRICT FROM BELOW WHEN SICN>0.15
!
         SICROW(I)=MAX(SICROW(I),SIC_MIN)
        ELSE
!
!       * DEFINE SIC WHEN SICN<=0.15
!
         SICROW(I)=(SIC_MIN/0.15)*SICNROW(I)
        ENDIF
       ENDIF
      ENDDO
!
      RETURN
    END SUBROUTINE PHSIL8

end module atm_f_data
