#ifdef AIX
@PROCESS ALIAS_SIZE(805306368)
#endif

module dice_comp_mod 1,16

! !USES:

  use shr_const_mod
  use shr_sys_mod
  use shr_kind_mod , only: IN=>SHR_KIND_IN, R8=>SHR_KIND_R8, &
                           CS=>SHR_KIND_CS, CL=>SHR_KIND_CL
  use shr_file_mod , only: shr_file_getunit, shr_file_getlogunit, shr_file_getloglevel, &
                           shr_file_setlogunit, shr_file_setloglevel, shr_file_setio, &
                           shr_file_freeunit
  use shr_mpi_mod  , only: shr_mpi_bcast
  use shr_flux_mod , only: shr_flux_atmIce
  use shr_cal_mod  , only: shr_cal_ymd2eday
  use mct_mod
  use esmf_mod
  use perf_mod

  use shr_strdata_mod
  use shr_dmodel_mod
  use shr_pcdf_mod

  use seq_cdata_mod
  use seq_infodata_mod
  use seq_timemgr_mod
  use seq_flds_indices , only: nflds_i2x, &
                               nflds_x2i
  use seq_flds_mod     , only: seq_flds_i2x_fields, &
                               seq_flds_x2i_fields
!
! !PUBLIC TYPES:
  implicit none
  private ! except

!--------------------------------------------------------------------------
! Public interfaces
!--------------------------------------------------------------------------

  public :: dice_comp_init
  public :: dice_comp_run
  public :: dice_comp_final

!--------------------------------------------------------------------------
! Private data
!--------------------------------------------------------------------------

  !--- other ---
  character(CS) :: myModelName = 'ice'   ! user defined model name
  integer(IN)   :: mpicom
  integer(IN)   :: my_task               ! my task in mpi communicator mpicom
  integer(IN)   :: npes                  ! total number of tasks
  integer(IN),parameter :: master_task=0 ! task number of master task
  integer(IN)   :: logunit               ! logging unit number
  character(CL) :: ice_mode              ! mode
  integer(IN)   :: dbug = 0              ! debug level (higher is more)
  logical       :: firstcall             ! first call logical
  logical       :: scmMode = .false.     ! single column mode
  real(R8)      :: scmLat  = shr_const_SPVAL  ! single column lat
  real(R8)      :: scmLon  = shr_const_SPVAL  ! single column lon
  logical       :: read_restart          ! start from restart
  real(R8)      :: flux_swpf             ! short-wave penatration factor
  real(R8)      :: flux_Qmin             ! bound on melt rate
  logical       :: flux_Qacc             ! activates water accumulation/melt wrt Q
  real(R8)      :: flux_Qacc0            ! initial water accumulation value

  character(len=*),parameter :: rpfile = 'rpointer.ice'
  character(len=*),parameter :: nullstr = 'undefined'

  real(R8),parameter  :: pi     = shr_const_pi      ! pi
  real(R8),parameter  :: spval  = shr_const_spval   ! flags invalid data
  real(R8),parameter  :: tFrz   = shr_const_tkfrzsw ! temp of freezing salt-water
  real(R8),parameter  :: latice = shr_const_latice  ! latent heat of fusion
  real(R8),parameter  :: cDay   = shr_const_cDay    ! sec in calendar day
  real(R8),parameter  :: waterMax = 1000.0_R8        ! wrt iFrac comp & frazil ice (kg/m^2)

  !----- surface albedo constants ------
  real(R8),parameter  :: snwfrac = 0.286_R8 ! snow cover fraction ~ [0,1]
  real(R8),parameter  :: as_nidf = 0.950_R8 ! albedo: snow,near-infr,diffuse
  real(R8),parameter  :: as_vsdf = 0.700_R8 ! albedo: snow,visible  ,diffuse
  real(R8),parameter  :: as_nidr = 0.960_R8 ! albedo: snow,near-infr,direct
  real(R8),parameter  :: as_vsdr = 0.800_R8 ! albedo: snow,visible  ,direct
  real(R8),parameter  :: ai_nidf = 0.700_R8 ! albedo: ice, near-infr,diffuse
  real(R8),parameter  :: ai_vsdf = 0.500_R8 ! albedo: ice, visible  ,diffuse
  real(R8),parameter  :: ai_nidr = 0.700_R8 ! albedo: ice, near-infr,direct
  real(R8),parameter  :: ai_vsdr = 0.500_R8 ! albedo: ice, visible  ,direct
  real(R8),parameter  :: ax_nidf = ai_nidf*(1.0_R8-snwfrac) + as_nidf*snwfrac
  real(R8),parameter  :: ax_vsdf = ai_vsdf*(1.0_R8-snwfrac) + as_vsdf*snwfrac
  real(R8),parameter  :: ax_nidr = ai_nidr*(1.0_R8-snwfrac) + as_nidr*snwfrac
  real(R8),parameter  :: ax_vsdr = ai_vsdr*(1.0_R8-snwfrac) + as_vsdr*snwfrac

  integer(IN) :: kswvdr,kswndr,kswvdf,kswndf,kq,kz,kua,kva,kptem,kshum,kdens,ktbot
  integer(IN) :: kiFrac,kt,kavsdr,kanidr,kavsdf,kanidf,kswnet,kmelth,kmeltw
  integer(IN) :: ksen,klat,klwup,kevap,ktauxa,ktauya,ktref,kqref,kswpen,ktauxo,ktauyo,ksalt

  type(shr_strdata_type) :: SDICE
  type(mct_rearr) :: rearr
!  type(mct_avect) :: avstrm   ! av of data from stream
  integer(IN) , pointer :: imask(:)
  real(R8)    , pointer :: yc(:)
  real(R8)    , pointer :: water(:)
!  real(R8)    , pointer :: ifrac0(:)

  integer(IN),parameter :: ktrans = 42
  character(16),parameter  :: avofld(1:ktrans) = &
     (/"So_t            ","So_s            ","So_u            ","So_v            ", &
       "So_dhdx         ","So_dhdy         ","Fioo_q          ","Sa_z            ", &
       "Sa_u            ","Sa_v            ","Sa_ptem         ","Sa_tbot         ", &
       "Sa_shum         ","Sa_dens         ","Faxa_swndr      ","Faxa_swvdr      ", &
       "Faxa_swndf      ","Faxa_swvdf      ","Faxa_lwdn       ","Faxa_rain       ", &
       "Faxa_snow       ","Si_t            ","Si_tref         ","Si_qref         ", &
       "Si_ifrac        ","Si_avsdr        ","Si_anidr        ","Si_avsdf        ", &
       "Si_anidf        ","Faii_taux       ","Faii_tauy       ","Faii_lat        ", &
       "Faii_sen        ","Faii_lwup       ","Faii_evap       ","Faii_swnet      ", &
       "Fioi_swpen      ","Fioi_melth      ","Fioi_meltw      ","Fioi_salt       ", &
       "Fioi_taux       ","Fioi_tauy       " /)

  character(16),parameter  :: avifld(1:ktrans) = &
     (/"to              ","s               ","uo              ","vo              ", &
       "dhdx            ","dhdy            ","q               ","z               ", &
       "ua              ","va              ","ptem            ","tbot            ", &
       "shum            ","dens            ","swndr           ","swvdr           ", &
       "swndf           ","swvdf           ","lwdn            ","rain            ", &
       "snow            ","t               ","tref            ","qref            ", &
       "ifrac           ","avsdr           ","anidr           ","avsdf           ", &
       "anidf           ","tauxa           ","tauya           ","lat             ", &
       "sen             ","lwup            ","evap            ","swnet           ", &
       "swpen           ","melth           ","meltw           ","salt            ", &
       "tauxo           ","tauyo           " /)

  save

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CONTAINS
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: dice_comp_init
!
! !DESCRIPTION:
!     initialize data ice model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------


subroutine dice_comp_init( EClock, cdata, x2i, i2x, NLFilename ) 2,35

    implicit none

! !INPUT/OUTPUT PARAMETERS:

    type(ESMF_Clock)            , intent(in)    :: EClock
    type(seq_cdata)             , intent(inout) :: cdata
    type(mct_aVect)             , intent(inout) :: x2i, i2x
    character(len=*), optional  , intent(in)    :: NLFilename ! Namelist filename

!EOP

    !--- local variables ---
    integer(IN)   :: n,k         ! generic counters
    integer(IN)   :: ierr        ! error code
    integer(IN)   :: COMPID      ! comp id
    integer(IN)   :: gsize       ! global size
    integer(IN)   :: lsize     ! local size
    integer(IN)   :: shrlogunit, shrloglev ! original log unit and level
    integer(IN)   :: nunit       ! unit number
    integer(IN)   :: kfld        ! field reference
    logical       :: ice_present    ! flag
    logical       :: ice_prognostic ! flag

    type(seq_infodata_type), pointer :: infodata
    type(mct_gsMap)        , pointer :: gsmap
    type(mct_gGrid)        , pointer :: ggrid

    character(CL) :: filePath    ! generic file path
    character(CL) :: fileName    ! generic file name
    character(CS) :: timeName    ! domain file: time variable name
    character(CS) ::  lonName    ! domain file: lon  variable name
    character(CS) ::  latName    ! domain file: lat  variable name
    character(CS) :: maskName    ! domain file: mask variable name
    character(CS) :: areaName    ! domain file: area variable name

    integer(IN)   :: yearFirst   ! first year to use in data stream
    integer(IN)   :: yearLast    ! last  year to use in data stream
    integer(IN)   :: yearAlign   ! data year that aligns with yearFirst

    character(CL) :: ice_in      ! dshr ice namelist
    character(CL) :: decomp      ! decomp strategy
    character(CL) :: rest_file   ! restart filename
    character(CL) :: rest_file_strm   ! restart filename for stream
    character(CL) :: restfilm    ! model restart file namelist
    character(CL) :: restfils    ! stream restart file namelist
    logical       :: exists      ! file existance logical
    integer(IN)   :: nu          ! unit number

    !----- define namelist -----
    namelist / dice_nml / &
        ice_in, decomp, flux_swpf, flux_Qmin, flux_Qacc, flux_Qacc0, restfilm, restfils

    !--- formats ---
    character(*), parameter :: F00   = "('(dice_comp_init) ',8a)"
    character(*), parameter :: F01   = "('(dice_comp_init) ',a,5i8)"
    character(*), parameter :: F02   = "('(dice_comp_init) ',a,4es13.6)"
    character(*), parameter :: F03   = "('(dice_comp_init) ',a,i8,a)"
    character(*), parameter :: F04   = "('(dice_comp_init) ',2a,2i8,'s')"
    character(*), parameter :: F05   = "('(dice_comp_init) ',a,2f10.4)"
    character(*), parameter :: F06   = "('(dice_comp_init) ',a,5l3)"
    character(*), parameter :: F90   = "('(dice_comp_init) ',73('='))"
    character(*), parameter :: F91   = "('(dice_comp_init) ',73('-'))"
    character(*), parameter :: subName = "(dice_comp_init) "
!-------------------------------------------------------------------------------


    call t_startf('DICE_INIT')

    firstcall = .true.

    ! Set cdata pointers

    call seq_cdata_setptrs(cdata, ID=COMPID, mpicom=mpicom, &
         gsMap=gsmap, dom=ggrid, infodata=infodata)

    ! Determine communicator groups and sizes

    call mpi_comm_rank(mpicom, my_task, ierr)
    call mpi_comm_size(mpicom, npes, ierr)

    !--- open log file ---
    if (my_task == master_task) then
       logUnit = shr_file_getUnit()
       call shr_file_setIO('ice_modelio.nml',logUnit)
    else
       logUnit = 6
    endif

    !----------------------------------------------------------------------------
    ! Reset shr logging to my log file
    !----------------------------------------------------------------------------
    call shr_file_getLogUnit (shrlogunit)
    call shr_file_getLogLevel(shrloglev)
    call shr_file_setLogUnit (logUnit)

    !----------------------------------------------------------------------------
    ! Set a Few Defaults
    !----------------------------------------------------------------------------

    call seq_infodata_getData(infodata,single_column=scmMode, &
   &                          scmlat=scmlat, scmlon=scmLon)

    ice_present = .false.
    ice_prognostic = .false.
    call seq_infodata_GetData(infodata,read_restart=read_restart)

    !----------------------------------------------------------------------------
    ! Read dice_in
    !----------------------------------------------------------------------------

    call t_startf('dice_readnml')

    filename = "dice_in"
    ice_in = "unset"
    decomp = "1d"
    flux_swpf  =     0.0_R8  ! no penetration
    flux_Qmin  =  -300.0_R8  ! kg/s/m^2
    flux_Qacc  = .false.     ! no accumulation 
    flux_Qacc0 =     0.0_R8  ! no water
    restfilm = trim(nullstr)
    restfils = trim(nullstr)
    if (my_task == master_task) then
       nunit = shr_file_getUnit() ! get unused unit number
       open (nunit,file=trim(filename),status="old",action="read")
       read (nunit,nml=dice_nml,iostat=ierr)
       close(nunit)
       call shr_file_freeUnit(nunit)
       if (ierr > 0) then
          write(logunit,F01) 'ERROR: reading input namelist, '//trim(filename)//' iostat=',ierr
          call shr_sys_abort(subName//': namelist read error '//trim(filename))
       end if
       write(logunit,F00)' ice_in     = ',trim(ice_in)
       write(logunit,F00)' decomp     = ',trim(decomp)
       write(logunit,F02)' flux_swpf  = ',flux_swpf
       write(logunit,F02)' flux_Qmin  = ',flux_Qmin
       write(logunit,F06)' flux_Qacc  = ',flux_Qacc
       write(logunit,F02)' flux_Qacc0 = ',flux_Qacc0
       write(logunit,F00)' restfilm   = ',trim(restfilm)
       write(logunit,F00)' restfils   = ',trim(restfils)
    endif
    call shr_mpi_bcast(ice_in    ,mpicom,'ice_in')
    call shr_mpi_bcast(decomp    ,mpicom,'decomp')
    call shr_mpi_bcast(flux_swpf ,mpicom,'flux_swpf')
    call shr_mpi_bcast(flux_Qmin ,mpicom,'flux_Qmin')
    call shr_mpi_bcast(flux_Qacc ,mpicom,'flux_Qacc')
    call shr_mpi_bcast(flux_Qacc0,mpicom,'flux_Qacc0')
    call shr_mpi_bcast(restfilm,mpicom,'restfilm')
    call shr_mpi_bcast(restfils,mpicom,'restfils')

    rest_file = trim(restfilm)
    rest_file_strm = trim(restfils)

    !----------------------------------------------------------------------------
    ! Read dshr namelist
    !----------------------------------------------------------------------------

    call shr_strdata_readnml(SDICE,trim(ice_in),mpicom=mpicom)

    !----------------------------------------------------------------------------
    ! Validate mode
    !----------------------------------------------------------------------------

    ice_mode = trim(SDICE%dataMode)

    ! check that we know how to handle the mode

    if (trim(ice_mode) == 'NULL' .or. &
        trim(ice_mode) == 'SSTDATA' .or. &
        trim(ice_mode) == 'COPYALL') then
      if (my_task == master_task) &
         write(logunit,F00) ' ice mode = ',trim(ice_mode)
    else
      write(logunit,F00) ' ERROR illegal ice mode = ',trim(ice_mode)
      call shr_sys_abort()
    endif

    call t_stopf('dice_readnml')

    !----------------------------------------------------------------------------
    ! Initialize datasets
    !----------------------------------------------------------------------------

    call t_startf('dice_strdata_init')

    if (trim(ice_mode) /= 'NULL') then
       ice_present = .true.
       if (scmmode) then
          if (my_task == master_task) &
             write(logunit,F05) ' scm lon lat = ',scmlon,scmlat
          call shr_strdata_init(SDICE,mpicom,compid,name='ice', &
                      scmmode=scmmode,scmlon=scmlon,scmlat=scmlat)
       else
          call shr_strdata_init(SDICE,mpicom,compid,name='ice')
       endif
    endif

    if (trim(ice_mode) == 'SSTDATA' .or. &
        trim(ice_mode) == 'COPYALL') then
       ice_prognostic = .true.
    endif

    if (my_task == master_task) then
       call shr_strdata_print(SDICE,'SDICE data')
    endif

    call t_stopf('dice_strdata_init')

    !----------------------------------------------------------------------------
    ! Set flag to specify data components
    !----------------------------------------------------------------------------

    call seq_infodata_PutData(infodata, &
      ice_present=ice_present, ice_prognostic=ice_prognostic, &
      ice_nx=SDICE%nxg, ice_ny=SDICE%nyg )

    !----------------------------------------------------------------------------
    ! Initialize MCT global seg map, 1d decomp
    !----------------------------------------------------------------------------

    call t_startf('dice_initgsmaps')
    if (my_task == master_task) write(logunit,F00) ' initialize gsmaps'
    call shr_sys_flush(logunit)

    call shr_dmodel_gsmapcreate(gsmap,SDICE%nxg*SDICE%nyg,compid,mpicom,decomp)
    lsize = mct_gsmap_lsize(gsmap,mpicom)

    if (ice_present) then
       call mct_rearr_init(SDICE%gsmap,gsmap,mpicom,rearr)
    endif

    call t_stopf('dice_initgsmaps')

    !----------------------------------------------------------------------------
    ! Initialize MCT domain
    !----------------------------------------------------------------------------

    call t_startf('dice_initmctdom')
    if (my_task == master_task) write(logunit,F00) 'copy domains'
    call shr_sys_flush(logunit)

    if (ice_present) call shr_dmodel_rearrGGrid(SDICE%grid, ggrid, gsmap, rearr, mpicom)

    call t_stopf('dice_initmctdom')

    !----------------------------------------------------------------------------
    ! Initialize MCT attribute vectors
    !----------------------------------------------------------------------------

    call t_startf('dice_initmctavs')
    if (my_task == master_task) write(logunit,F00) 'allocate AVs'
    call shr_sys_flush(logunit)

    call mct_aVect_init(i2x, rList=seq_flds_i2x_fields, lsize=lsize)
    call mct_aVect_zero(i2x)

    kiFrac = mct_aVect_indexRA(i2x,'Si_ifrac')
    kt     = mct_aVect_indexRA(i2x,'Si_t')
    ktref  = mct_aVect_indexRA(i2x,'Si_tref')
    kqref  = mct_aVect_indexRA(i2x,'Si_qref')
    kavsdr = mct_aVect_indexRA(i2x,'Si_avsdr')
    kanidr = mct_aVect_indexRA(i2x,'Si_anidr')
    kavsdf = mct_aVect_indexRA(i2x,'Si_avsdf')
    kanidf = mct_aVect_indexRA(i2x,'Si_anidf')
    kswnet = mct_aVect_indexRA(i2x,'Faii_swnet')
    ksen   = mct_aVect_indexRA(i2x,'Faii_sen')
    klat   = mct_aVect_indexRA(i2x,'Faii_lat')
    klwup  = mct_aVect_indexRA(i2x,'Faii_lwup')
    kevap  = mct_aVect_indexRA(i2x,'Faii_evap')
    ktauxa = mct_aVect_indexRA(i2x,'Faii_taux')
    ktauya = mct_aVect_indexRA(i2x,'Faii_tauy')
    kmelth = mct_aVect_indexRA(i2x,'Fioi_melth')
    kmeltw = mct_aVect_indexRA(i2x,'Fioi_meltw')
    kswpen = mct_aVect_indexRA(i2x,'Fioi_swpen')
    ktauxo = mct_aVect_indexRA(i2x,'Fioi_taux')
    ktauyo = mct_aVect_indexRA(i2x,'Fioi_tauy')
    ksalt  = mct_aVect_indexRA(i2x,'Fioi_salt')

    call mct_aVect_init(x2i, rList=seq_flds_x2i_fields, lsize=lsize)
    call mct_aVect_zero(x2i)

    kswvdr = mct_aVect_indexRA(x2i,'Faxa_swvdr')
    kswndr = mct_aVect_indexRA(x2i,'Faxa_swndr')
    kswvdf = mct_aVect_indexRA(x2i,'Faxa_swvdf')
    kswndf = mct_aVect_indexRA(x2i,'Faxa_swndf')
    kq     = mct_aVect_indexRA(x2i,'Fioo_q')
    kz     = mct_aVect_indexRA(x2i,'Sa_z')
    kua    = mct_aVect_indexRA(x2i,'Sa_u')
    kva    = mct_aVect_indexRA(x2i,'Sa_v')
    kptem  = mct_aVect_indexRA(x2i,'Sa_ptem')
    kshum  = mct_aVect_indexRA(x2i,'Sa_shum')
    kdens  = mct_aVect_indexRA(x2i,'Sa_dens')
    ktbot  = mct_aVect_indexRA(x2i,'Sa_tbot')

!    call mct_aVect_init(avstrm, rList=flds_strm, lsize=lsize)
!    call mct_aVect_zero(avstrm)

    allocate(imask(lsize))
    allocate(yc(lsize))
    allocate(water(lsize))
!    allocate(iFrac0(lsize))

    kfld = mct_aVect_indexRA(ggrid%data,'mask')
    imask(:) = nint(ggrid%data%rAttr(kfld,:))
    kfld = mct_aVect_indexRA(ggrid%data,'lat')
    yc(:) = ggrid%data%rAttr(kfld,:)

    call t_stopf('dice_initmctavs')

    !----------------------------------------------------------------------------
    ! Read restart
    !----------------------------------------------------------------------------

    if (read_restart) then
       if (trim(rest_file)      == trim(nullstr) .and. &
           trim(rest_file_strm) == trim(nullstr)) then
          if (my_task == master_task) then
             write(logunit,F00) ' restart filenames from rpointer'
             call shr_sys_flush(logunit)
             inquire(file=trim(rpfile),exist=exists)
             if (.not.exists) then
                write(logunit,F00) ' ERROR: rpointer file does not exist'
                call shr_sys_abort(trim(subname)//' ERROR: rpointer file missing')
             endif
             nu = shr_file_getUnit()
             open(nu,file=trim(rpfile),form='formatted')
             read(nu,'(a)') rest_file
             read(nu,'(a)') rest_file_strm
             close(nu)
             call shr_file_freeUnit(nu)
             inquire(file=trim(rest_file_strm),exist=exists)
          endif
          call shr_mpi_bcast(rest_file,mpicom,'rest_file')
          call shr_mpi_bcast(rest_file_strm,mpicom,'rest_file_strm')
       else
          ! use namelist already read
          if (my_task == master_task) then
             write(logunit,F00) ' restart filenames from namelist '
             call shr_sys_flush(logunit)
             inquire(file=trim(rest_file_strm),exist=exists)
          endif
       endif
       call shr_mpi_bcast(exists,mpicom,'exists')
       if (my_task == master_task) write(logunit,F00) ' reading ',trim(rest_file)
       call shr_pcdf_readwrite('read',trim(rest_file),mpicom,gsmap,rf1=water,rf1n='water')
       if (exists) then
          if (my_task == master_task) write(logunit,F00) ' reading ',trim(rest_file_strm)
          call shr_strdata_restRead(trim(rest_file_strm),SDICE,mpicom)
       else
          if (my_task == master_task) write(logunit,F00) ' file not found, skipping ',trim(rest_file_strm)
       endif
       call shr_sys_flush(logunit)
    endif

    !----------------------------------------------------------------------------
    ! On initial call, x2i is unset, so set for use in run method
    !  These values should have no impact on the solution!!
    !----------------------------------------------------------------------------
    x2i%rAttr(kz,:)    = 10.0_R8
    x2i%rAttr(kua,:)   = 5.0_R8
    x2i%rAttr(kva,:)   = 5.0_R8
    x2i%rAttr(kptem,:) = 260.0_R8
    x2i%rAttr(ktbot,:) = 260.0_R8
    x2i%rAttr(kshum,:) = 0.0014_R8
    x2i%rAttr(kdens,:) = 1.3_R8

    !----------------------------------------------------------------------------
    ! Set initial ice state, needed for CCSM atm initialization
    !----------------------------------------------------------------------------

    call t_adj_detailf(+2)
    call dice_comp_run( EClock, cdata,  x2i, i2x)
    call t_adj_detailf(-2)

    !----------------------------------------------------------------------------
    ! Reset shr logging to original values
    !----------------------------------------------------------------------------

    if (my_task == master_task) write(logunit,F00) 'dice_comp_init done'
    call shr_sys_flush(logunit)

    call shr_file_setLogUnit (shrlogunit)
    call shr_file_setLogLevel(shrloglev)
    call shr_sys_flush(logunit)

    call t_stopf('DICE_INIT')

end subroutine dice_comp_init

!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: dice_comp_run
!
! !DESCRIPTION:
!     run method for dead ice model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------


subroutine dice_comp_run( EClock, cdata,  x2i, i2x) 2,24

   implicit none

! !INPUT/OUTPUT PARAMETERS:

   type(ESMF_Clock)            ,intent(in)    :: EClock
   type(seq_cdata)             ,intent(inout) :: cdata
   type(mct_aVect)             ,intent(inout) :: x2i        ! driver -> dead
   type(mct_aVect)             ,intent(inout) :: i2x        ! dead   -> driver

!EOP

   !--- local ---
   type(mct_gsMap)        , pointer :: gsmap
   type(mct_gGrid)        , pointer :: ggrid

   integer(IN)   :: CurrentYMD        ! model date
   integer(IN)   :: CurrentTOD        ! model sec into model date
   integer(IN)   :: yy,mm,dd          ! year month day
   integer(IN)   :: n                 ! indices
   integer(IN)   :: nf                ! fields loop index
   integer(IN)   :: nl                ! ice frac index
   integer(IN)   :: lsize           ! size of attr vect
   integer(IN)   :: shrlogunit, shrloglev ! original log unit and level
   logical       :: glcrun_alarm      ! is glc going to run now
   logical       :: newdata           ! has newdata been read
   logical       :: mssrmlf           ! remove old data
   integer(IN)   :: idt               ! integer timestep
   real(R8)      :: dt                ! timestep
   real(R8)      :: hn                ! h field
   logical       :: write_restart     ! restart now
   character(CL) :: case_name         ! case name
   character(CL) :: rest_file         ! restart_file
   character(CL) :: rest_file_strm    ! restart_file for stream
   integer(IN)   :: nu                ! unit number
   real(R8)      :: qmeltall          ! q that would melt all accumulated water
   real(R8)      :: cosarg            ! for setting ice temp pattern
   integer(IN)   :: eday, eday0       ! elapsed day counters

   type(seq_infodata_type), pointer :: infodata

   character(*), parameter :: F00   = "('(dice_comp_run) ',8a)"
   character(*), parameter :: F04   = "('(dice_comp_run) ',2a,2i8,'s')"
   character(*), parameter :: subName = "(dice_comp_run) "
!-------------------------------------------------------------------------------

   call t_startf('DICE_RUN')

   call t_startf('dice_run1')

  !----------------------------------------------------------------------------
  ! Reset shr logging to my log file
  !----------------------------------------------------------------------------
   call shr_file_getLogUnit (shrlogunit)
   call shr_file_getLogLevel(shrloglev)
   call shr_file_setLogUnit (logUnit)

   call seq_cdata_setptrs(cdata, gsMap=gsmap, dom=ggrid)

   call seq_cdata_setptrs(cdata, infodata=infodata)

   call seq_timemgr_EClockGetData( EClock, curr_ymd=CurrentYMD, curr_tod=CurrentTOD)
   call seq_timemgr_EClockGetData( EClock, curr_yr=yy, curr_mon=mm, curr_day=dd)
   call seq_timemgr_EClockGetData( EClock, dtime=idt)
   dt = idt * 1.0_r8
   write_restart = seq_timemgr_RestartAlarmIsOn(EClock)

   call t_stopf('dice_run1')

   !--------------------
   ! UNPACK
   !--------------------

   call t_startf('dice_unpack')

!  lsize = mct_avect_lsize(x2i)

!   do nf=1,nflds_x2i
!   do n=1,lsize
!     ?? = x2i%rAttr(nf,n)
!   enddo
!   enddo

   call t_stopf('dice_unpack')

   !--------------------
   ! ADVANCE ICE
   !--------------------

   call t_barrierf('dice_BARRIER',mpicom)
   call t_startf('dice')

   !--- copy all fields from streams to i2x as default ---

   if (trim(ice_mode) /= 'NULL') then
      call t_startf('dice_strdata_advance')
      call shr_strdata_advance(SDICE,currentYMD,currentTOD,mpicom,'dice')
      call t_stopf('dice_strdata_advance')
      call t_barrierf('dice_scatter_BARRIER',mpicom)
      call t_startf('dice_scatter')
      do n = 1,SDICE%nstreams
         call shr_dmodel_translateAV(SDICE%avs(n),i2x,avifld,avofld,rearr)
      enddo
      call t_stopf('dice_scatter')
   else
      call mct_aVect_zero(i2x)
   endif

   call t_startf('dice_mode')

   select case (trim(ice_mode))

   case('COPYALL') 
      ! do nothing extra

   case('SSTDATA')
      if (firstcall .and. .not. read_restart) then
!         iFrac0 = iFrac  ! previous step's ice fraction
         water  = 0.0_R8 ! previous step's water accumulation
         where (i2x%rAttr(kiFrac,:) > 0.0_R8) water(:) = flux_Qacc0 
      endif

      call shr_cal_ymd2eday(0,mm,dd,eDay )    ! model date
      call shr_cal_ymd2eday(0,09,01,eDay0)    ! sept 1st
      cosArg = 2.0_R8*pi*(real(eDay,R8) + real(currentTOD,R8)/cDay - real(eDay0,R8))/365.0_R8

      lsize = mct_avect_lsize(i2x)

      do n = 1,lsize

         !--- fix erroneous iFrac ---
         i2x%rAttr(kiFrac,n) = min(1.0_R8,max(0.0_R8,i2x%rAttr(kiFrac,n))) 

         !--- fabricate ice surface T, fix erroneous iFrac ---
         if ( yc(n) > 0.0_R8) then 
            i2x%rAttr(kt,n) = 260.0_R8 + 10.0_R8*cos(cosArg)
         else
            i2x%rAttr(kt,n) = 260.0_R8 - 10.0_R8*cos(cosArg)
         end if

         !--- set albedos (constant) ---
         i2x%rAttr(kavsdr,n) = ax_vsdr 
         i2x%rAttr(kanidr,n) = ax_nidr 
         i2x%rAttr(kavsdf,n) = ax_vsdf 
         i2x%rAttr(kanidf,n) = ax_nidf 

         !--- swnet is sent to cpl as a diagnostic quantity only ---
         !--- newly recv'd swdn goes with previously sent albedo ---
         !--- but albedos are (currently) time invariant         ---
         i2x%rAttr(kswnet,n) = (1.0_R8 - i2x%rAttr(kavsdr,n))*x2i%rAttr(kswvdr,n) &
         &                   + (1.0_R8 - i2x%rAttr(kanidr,n))*x2i%rAttr(kswndr,n) &
         &                   + (1.0_R8 - i2x%rAttr(kavsdf,n))*x2i%rAttr(kswvdf,n) &
         &                   + (1.0_R8 - i2x%rAttr(kanidf,n))*x2i%rAttr(kswndf,n)

         !--- compute melt/freeze water balance, adjust iFrac  -------------
         if ( .not. flux_Qacc ) then ! Q accumulation option is OFF
            i2x%rAttr(kmelth,n) = min(x2i%rAttr(kq,n),0.0_R8 ) ! q<0 => melt potential
            i2x%rAttr(kmelth,n) = max(i2x%rAttr(kmelth,n),Flux_Qmin   ) ! limit the melt rate
            i2x%rAttr(kmeltw,n) =    -i2x%rAttr(kmelth,n)/latice   ! corresponding water flux

         else                                 ! Q accumulation option is ON
            !--------------------------------------------------------------
            ! 1a) Q<0 & iFrac > 0  =>  infinite supply of water to melt
            ! 1b) Q<0 & iFrac = 0  =>  melt accumulated water only
            ! 2a) Q>0 & iFrac > 0  =>  zero-out accumulated water
            ! 2b) Q>0 & iFrac = 0  =>  accumulated water
            !--------------------------------------------------------------
            if ( x2i%rAttr(kq,n) <  0.0_R8 ) then ! Q<0 => melt
               if (i2x%rAttr(kiFrac,n) > 0.0_R8 ) then
                  i2x%rAttr(kmelth,n) = i2x%rAttr(kiFrac,n)*max(x2i%rAttr(kq,n),Flux_Qmin)
                  i2x%rAttr(kmeltw,n) =    -i2x%rAttr(kmelth,n)/latice
               !  water(n) = < don't change this value >
               else
                  Qmeltall   = -water(n)*latice/dt
                  i2x%rAttr(kmelth,n) = max(x2i%rAttr(kq,n), Qmeltall, Flux_Qmin )
                  i2x%rAttr(kmeltw,n) = -i2x%rAttr(kmelth,n)/latice
                  water(n) =  water(n) - i2x%rAttr(kmeltw,n)*dt
               end if
            else                       ! Q>0 => freeze
               if (i2x%rAttr(kiFrac,n) > 0.0_R8 ) then
                  i2x%rAttr(kmelth,n) = 0.0_R8
                  i2x%rAttr(kmeltw,n) = 0.0_R8
                  water(n) = 0.0_R8
               else
                  i2x%rAttr(kmelth,n) = 0.0_R8
                  i2x%rAttr(kmeltw,n) = 0.0_R8
                  water(n) = water(n) + dt*i2x%rAttr(kq,n)/latice
               end if
            end if

            if (water(n) < 1.0e-16_R8 ) water(n) = 0.0_R8

            !--- non-zero water => non-zero iFrac ---
            if (i2x%rAttr(kiFrac,n) <= 0.0_R8  .and.  water(n) > 0.0_R8) then
               i2x%rAttr(kiFrac,n) = min(1.0_R8,water(n)/waterMax)
               ! i2x%rAttr(kT,n) = Tfrz     ! T can be above freezing?!?
            end if

            !--- cpl multiplies melth & meltw by iFrac ---
            !--- divide by iFrac here => fixed quantity flux (not per area) ---
            if (i2x%rAttr(kiFrac,n) > 0.0_R8) then
               i2x%rAttr(kiFrac,n) = max( 0.01_R8, i2x%rAttr(kiFrac,n)) ! min iFrac
               i2x%rAttr(kmelth,n) = i2x%rAttr(kmelth,n)/i2x%rAttr(kiFrac,n)
               i2x%rAttr(kmeltw,n) = i2x%rAttr(kmeltw,n)/i2x%rAttr(kiFrac,n)
            else
               i2x%rAttr(kmelth,n) = 0.0_R8
               i2x%rAttr(kmeltw,n) = 0.0_R8
            end if
         end if

         !--- modify T wrt iFrac: (iFrac -> 0) => (T -> Tfrz) ---
         i2x%rAttr(kt,n) = Tfrz + i2x%rAttr(kiFrac,n)*(i2x%rAttr(kt,n)-Tfrz) 

      end do

      !----------------------------------------------------------------------------
      ! compute atm/ice surface fluxes
      !----------------------------------------------------------------------------
      call shr_flux_atmIce(iMask  ,x2i%rAttr(kz,:)     ,x2i%rAttr(kua,:)    ,x2i%rAttr(kva,:), &
               x2i%rAttr(kptem,:) ,x2i%rAttr(kshum,:)  ,x2i%rAttr(kdens,:)  ,x2i%rAttr(ktbot,:),  &
               i2x%rAttr(kt,:)    ,i2x%rAttr(ksen,:)   ,i2x%rAttr(klat,:)   ,i2x%rAttr(klwup,:), &
               i2x%rAttr(kevap,:) ,i2x%rAttr(ktauxa,:) ,i2x%rAttr(ktauya,:) ,i2x%rAttr(ktref,:), &
               i2x%rAttr(kqref,:) )

      !----------------------------------------------------------------------------
      ! compute ice/oce surface fluxes (except melth & meltw, see above)
      !----------------------------------------------------------------------------
      do n=1,lsize
         if (iMask(n) == 0) then
            i2x%rAttr(kswpen,n) = spval
            i2x%rAttr(kmelth,n) = spval
            i2x%rAttr(kmeltw,n) = spval
            i2x%rAttr(ksalt ,n) = spval
            i2x%rAttr(ktauxo,n) = spval
            i2x%rAttr(ktauyo,n) = spval
            i2x%rAttr(kiFrac,n) = 0.0_R8
         else
            !--- penetrating short wave ---
            i2x%rAttr(kswpen,n) = max(0.0_R8, flux_swpf*i2x%rAttr(kswnet,n) ) ! must be non-negative

            !--- i/o surface stress ( = atm/ice stress) ---
            i2x%rAttr(ktauxo,n) = i2x%rAttr(ktauxa,n)
            i2x%rAttr(ktauyo,n) = i2x%rAttr(ktauya,n)

            !--- salt flux ---
           i2x%rAttr(ksalt ,n) = 0.0_R8
         end if

!         !--- save ifrac for next timestep
!         iFrac0(n) = i2x%rAttr(kiFrac,n)
      end do


   end select

   call t_stopf('dice_mode')

   if (write_restart) then
      call t_startf('dice_restart')
      call seq_infodata_GetData( infodata, case_name=case_name)
      write(rest_file,"(2a,i4.4,a,i2.2,a,i2.2,a,i5.5,a)") &
        trim(case_name), '.dice.r.', yy,'-',mm,'-',dd,'-',currentTOD,'.nc'
      write(rest_file_strm,"(2a,i4.4,a,i2.2,a,i2.2,a,i5.5,a)") &
        trim(case_name), '.dice.rs1.', yy,'-',mm,'-',dd,'-',currentTOD,'.bin'
      if (my_task == master_task) then
         nu = shr_file_getUnit()
         open(nu,file=trim(rpfile),form='formatted')
         write(nu,'(a)') rest_file
         write(nu,'(a)') rest_file_strm
         close(nu)
         call shr_file_freeUnit(nu)
      endif
      if (my_task == master_task) write(logunit,F04) ' writing ',trim(rest_file),currentYMD,currentTOD
       call shr_pcdf_readwrite('write',trim(rest_file),mpicom,gsmap,clobber=.true., &
         rf1=water,rf1n='water')
      if (my_task == master_task) write(logunit,F04) ' writing ',trim(rest_file_strm),currentYMD,currentTOD
      call shr_strdata_restWrite(trim(rest_file_strm),SDICE,mpicom,trim(case_name),'SDICE strdata')
      call shr_sys_flush(logunit)
      call t_stopf('dice_restart')
   endif

   call t_stopf('dice')

   !----------------------------------------------------------------------------
   ! Log output for model date
   ! Reset shr logging to original values
   !----------------------------------------------------------------------------

   call t_startf('dice_run2')
   if (my_task == master_task) then
      write(logunit,F04) trim(myModelName),': model date ', CurrentYMD,CurrentTOD
      call shr_sys_flush(logunit)
   end if
   firstcall = .false.
      
   call shr_file_setLogUnit (shrlogunit)
   call shr_file_setLogLevel(shrloglev)
   call shr_sys_flush(logunit)
   call t_stopf('dice_run2')

   call t_stopf('DICE_RUN')

end subroutine dice_comp_run

!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: dice_comp_final
!
! !DESCRIPTION:
!     finalize method for dead ice model
!
! !REVISION HISTORY:
!
! !INTERFACE: ------------------------------------------------------------------
!

subroutine dice_comp_final() 1

   implicit none

!EOP

   !--- formats ---
   character(*), parameter :: F00   = "('(dice_comp_final) ',8a)"
   character(*), parameter :: F91   = "('(dice_comp_final) ',73('-'))"
   character(*), parameter :: subName = "(dice_comp_final) "

!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------

   call t_startf('DICE_FINAL')

   if (my_task == master_task) then
      write(logunit,F91) 
      write(logunit,F00) trim(myModelName),': end of main integration loop'
      write(logunit,F91) 
   end if
      
   call t_stopf('DICE_FINAL')

end subroutine dice_comp_final
!===============================================================================
!===============================================================================

end module dice_comp_mod