!  SVN:$Id: ice_restart_driver.F90 1304 2015-11-12 19:28:14Z afrobert@nps.edu $
!=======================================================================

! Read and write ice model restart files
!
! authors Elizabeth C. Hunke, LANL
!         William H. Lipscomb LANL
!         David Bailey, NCAR
!
! 2004-05: Block structure added by William Lipscomb
!          Restart module separated from history module
! 2006 ECH: Accepted some CCSM code into mainstream CICE
!           Converted to free source form (F90) 
! 2008 ECH: Rearranged order in which internal stresses are written and read
! 2010 ECH: Changed eice, esno to qice, qsno
! 2012 ECH: Added routines for reading/writing extended grid
! 2013 DAB: Added generic interfaces for writing restart fields.

      module ice_restart_driver

      use ice_kinds_mod
      use ice_restart_shared, only: &
          restart, restart_ext, restart_dir, restart_file, pointer_file, &
          runid, runtype, use_restart_time, restart_format, lcdf64, lenstr
      use ice_restart

      implicit none
      private
      public :: dumpfile, restartfile, restartfile_v4
      save

!=======================================================================

      contains

!=======================================================================

!=======================================================================
!---! these subroutines write/read Fortran unformatted data files ..
!=======================================================================

! Dumps all values needed for a restart
! author Elizabeth C. Hunke, LANL

      subroutine dumpfile(filename_spec)

      use ice_blocks, only: nx_block, ny_block
      use ice_calendar, only: sec, month, mday, nyr, istep1, &
                              time, time_forc, year_init
      use ice_communicate, only: my_task, master_task
      use ice_constants, only: c0, c1
      use ice_domain, only: nblocks
      use ice_domain_size, only: nilyr, nslyr, ncat, max_blocks
      use ice_fileunits, only: nu_diag, nu_rst_pointer, nu_dump
      use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, &
          strocnxT, strocnyT, sst, frzmlt, iceumask, coszen, &
          stressp_1, stressp_2, stressp_3, stressp_4, &
          stressm_1, stressm_2, stressm_3, stressm_4, &
          stress12_1, stress12_2, stress12_3, stress12_4
      use ice_ocean, only: oceanmixed_ice
      use ice_read_write, only: ice_open, ice_write
      use ice_state, only: aicen, vicen, vsnon, trcrn, &
          nt_Tsfc, nt_sice, nt_qice, nt_qsno, uvel, vvel

      character(len=char_len_long), intent(in), optional :: filename_spec

      ! local variables

      integer (kind=int_kind) :: &
          i, j, k, n, iblk, &     ! counting indices
          iyear, imonth, iday     ! year, month, day

      character(len=char_len_long) :: filename

      logical (kind=log_kind) :: diag

      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
         work1

      character (len=3) :: nchar

      if (present(filename_spec)) then
         call init_restart_write(filename_spec)
      else
         call init_restart_write
      endif

      diag = .true.

      !-----------------------------------------------------------------
      ! state variables
      ! Tsfc is the only tracer written to binary files.  All other
      ! tracers are written to their own dump/restart binary files.
      !-----------------------------------------------------------------

      call write_restart_field(nu_dump,0,aicen(:,:,:,:),'ruf8','aicen',ncat,diag)
      call write_restart_field(nu_dump,0,vicen(:,:,:,:),'ruf8','vicen',ncat,diag)
      call write_restart_field(nu_dump,0,vsnon(:,:,:,:),'ruf8','vsnon',ncat,diag)
      call write_restart_field(nu_dump,0,trcrn(:,:,nt_Tsfc,:,:),'ruf8','Tsfcn',ncat,diag)

      do k=1,nilyr
         write(nchar,'(i3.3)') k
         call write_restart_field(nu_dump,0,trcrn(:,:,nt_sice+k-1,:,:),'ruf8', &
                                 'sice'//trim(nchar),ncat,diag)
      enddo

      do k=1,nilyr
         write(nchar,'(i3.3)') k
         call write_restart_field(nu_dump,0,trcrn(:,:,nt_qice+k-1,:,:),'ruf8', &
                                 'qice'//trim(nchar),ncat,diag)
      enddo

      do k=1,nslyr
         write(nchar,'(i3.3)') k
         call write_restart_field(nu_dump,0,trcrn(:,:,nt_qsno+k-1,:,:),'ruf8', &
                                 'qsno'//trim(nchar),ncat,diag)
      enddo

      !-----------------------------------------------------------------
      ! velocity
      !-----------------------------------------------------------------
      call write_restart_field(nu_dump,0,uvel,'ruf8','uvel',1,diag)
      call write_restart_field(nu_dump,0,vvel,'ruf8','vvel',1,diag)

      !-----------------------------------------------------------------
      ! radiation fields
      !-----------------------------------------------------------------
#ifdef CCSMCOUPLED
      call write_restart_field(nu_dump,0,coszen,'ruf8','coszen',1,diag)
#endif
      call write_restart_field(nu_dump,0,scale_factor,'ruf8','scale_factor',1,diag)

      call write_restart_field(nu_dump,0,swvdr,'ruf8','swvdr',1,diag)
      call write_restart_field(nu_dump,0,swvdf,'ruf8','swvdf',1,diag)
      call write_restart_field(nu_dump,0,swidr,'ruf8','swidr',1,diag)
      call write_restart_field(nu_dump,0,swidf,'ruf8','swidf',1,diag)

      !-----------------------------------------------------------------
      ! ocean stress (for bottom heat flux in thermo)
      !-----------------------------------------------------------------
      call write_restart_field(nu_dump,0,strocnxT,'ruf8','strocnxT',1,diag)
      call write_restart_field(nu_dump,0,strocnyT,'ruf8','strocnyT',1,diag)

      !-----------------------------------------------------------------
      ! internal stress
      !-----------------------------------------------------------------
      call write_restart_field(nu_dump,0,stressp_1,'ruf8','stressp_1',1,diag)
      call write_restart_field(nu_dump,0,stressp_3,'ruf8','stressp_3',1,diag)
      call write_restart_field(nu_dump,0,stressp_2,'ruf8','stressp_2',1,diag)
      call write_restart_field(nu_dump,0,stressp_4,'ruf8','stressp_4',1,diag)

      call write_restart_field(nu_dump,0,stressm_1,'ruf8','stressm_1',1,diag)
      call write_restart_field(nu_dump,0,stressm_3,'ruf8','stressm_3',1,diag)
      call write_restart_field(nu_dump,0,stressm_2,'ruf8','stressm_2',1,diag)
      call write_restart_field(nu_dump,0,stressm_4,'ruf8','stressm_4',1,diag)

      call write_restart_field(nu_dump,0,stress12_1,'ruf8','stress12_1',1,diag)
      call write_restart_field(nu_dump,0,stress12_3,'ruf8','stress12_3',1,diag)
      call write_restart_field(nu_dump,0,stress12_2,'ruf8','stress12_2',1,diag)
      call write_restart_field(nu_dump,0,stress12_4,'ruf8','stress12_4',1,diag)

      !-----------------------------------------------------------------
      ! ice mask for dynamics
      !-----------------------------------------------------------------
      
      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
      do iblk = 1, nblocks
         do j = 1, ny_block
         do i = 1, nx_block
            work1(i,j,iblk) = c0
            if (iceumask(i,j,iblk)) work1(i,j,iblk) = c1
         enddo
         enddo
      enddo
      !$OMP END PARALLEL DO
      call write_restart_field(nu_dump,0,work1,'ruf8','iceumask',1,diag)

      ! for mixed layer model
      if (oceanmixed_ice) then
         call write_restart_field(nu_dump,0,sst,'ruf8','sst',1,diag)
         call write_restart_field(nu_dump,0,frzmlt,'ruf8','frzmlt',1,diag)
      endif

      end subroutine dumpfile

!=======================================================================

! Restarts from a dump
! author Elizabeth C. Hunke, LANL

      subroutine restartfile (ice_ic)

      use ice_boundary, only: ice_HaloUpdate_stress
      use ice_broadcast, only: broadcast_scalar
      use ice_blocks, only: nghost, nx_block, ny_block
      use ice_calendar, only: istep0, istep1, time, time_forc, calendar, npt
      use ice_communicate, only: my_task, master_task
      use ice_constants, only: c0, p5, &
          field_loc_center, field_loc_NEcorner, &
          field_type_scalar, field_type_vector
      use ice_domain, only: nblocks, distrb_info, halo_info
      use ice_domain_size, only: nilyr, nslyr, ncat, nx_global, ny_global, &
          max_ntrcr, max_blocks
      use ice_fileunits, only: nu_diag, nu_rst_pointer, nu_restart
      use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, &
          strocnxT, strocnyT, sst, frzmlt, iceumask, coszen, &
          stressp_1, stressp_2, stressp_3, stressp_4, &
          stressm_1, stressm_2, stressm_3, stressm_4, &
          stress12_1, stress12_2, stress12_3, stress12_4
      use ice_gather_scatter, only: scatter_global_stress
      use ice_grid, only: tmask, grid_type
      use ice_itd, only: aggregate
      use ice_ocean, only: oceanmixed_ice
      use ice_read_write, only: ice_open, ice_read, ice_read_global
      use ice_state, only: trcr_depend, aice, vice, vsno, trcr, &
          aice0, aicen, vicen, vsnon, trcrn, aice_init, &
          nt_Tsfc, nt_sice, nt_qice, nt_qsno, uvel, vvel

      character (*), optional :: ice_ic

      ! local variables

      integer (kind=int_kind) :: &
         i, j, k, n, iblk, &     ! counting indices
         iignore                 ! dummy variable

      real (kind=real_kind) :: &
         rignore                 ! dummy variable

      character(len=char_len_long) :: &
         filename, filename0

      logical (kind=log_kind) :: &
         diag

      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
         work1

      real (kind=dbl_kind), dimension(:,:), allocatable :: &
         work_g1, work_g2

      character (len=3) :: nchar

      call init_restart_read(ice_ic)

      diag = .true.

      !-----------------------------------------------------------------
      ! state variables
      ! Tsfc is the only tracer read in this file.  All other
      ! tracers are in their own dump/restart files.
      !-----------------------------------------------------------------
      if (my_task == master_task) &
           write(nu_diag,*) ' min/max area, vol ice, vol snow, Tsfc'

      call read_restart_field(nu_restart,0,aicen,'ruf8', &
              'aicen',ncat,diag,field_loc_center, field_type_scalar)
      call read_restart_field(nu_restart,0,vicen,'ruf8', &
              'vicen',ncat,diag,field_loc_center, field_type_scalar)
      call read_restart_field(nu_restart,0,vsnon,'ruf8', &
              'vsnon',ncat,diag,field_loc_center, field_type_scalar)
      call read_restart_field(nu_restart,0,trcrn(:,:,nt_Tsfc,:,:),'ruf8', &
              'Tsfcn',ncat,diag,field_loc_center, field_type_scalar)

      if (my_task == master_task) &
         write(nu_diag,*) 'min/max sice for each layer'
      do k=1,nilyr
         write(nchar,'(i3.3)') k
         call read_restart_field(nu_restart,0,trcrn(:,:,nt_sice+k-1,:,:),'ruf8', &
              'sice'//trim(nchar),ncat,diag,field_loc_center,field_type_scalar)
      enddo

      if (my_task == master_task) &
         write(nu_diag,*) 'min/max qice for each layer'
      do k=1,nilyr
         write(nchar,'(i3.3)') k
         call read_restart_field(nu_restart,0,trcrn(:,:,nt_qice+k-1,:,:),'ruf8', &
              'qice'//trim(nchar),ncat,diag,field_loc_center,field_type_scalar)
      enddo

      if (my_task == master_task) &
         write(nu_diag,*) 'min/max qsno for each layer'
      do k=1,nslyr
         write(nchar,'(i3.3)') k
         call read_restart_field(nu_restart,0,trcrn(:,:,nt_qsno+k-1,:,:),'ruf8', &
              'qsno'//trim(nchar),ncat,diag,field_loc_center,field_type_scalar)
      enddo

      !-----------------------------------------------------------------
      ! velocity
      !-----------------------------------------------------------------
      if (my_task == master_task) &
           write(nu_diag,*) 'min/max velocity components'

      call read_restart_field(nu_restart,0,uvel,'ruf8', &
           'uvel',1,diag,field_loc_NEcorner, field_type_vector)
      call read_restart_field(nu_restart,0,vvel,'ruf8', &
           'vvel',1,diag,field_loc_NEcorner, field_type_vector)

      !-----------------------------------------------------------------
      ! radiation fields
      !-----------------------------------------------------------------

      if (my_task == master_task) &
         write(nu_diag,*) 'radiation fields'

#ifdef CCSMCOUPLED
      call read_restart_field(nu_restart,0,coszen,'ruf8', &
           'coszen',1,diag, field_loc_center, field_type_scalar)
#endif
      call read_restart_field(nu_restart,0,scale_factor,'ruf8', &
           'scale_factor',1,diag, field_loc_center, field_type_scalar)
      call read_restart_field(nu_restart,0,swvdr,'ruf8', &
           'swvdr',1,diag,field_loc_center, field_type_scalar)
      call read_restart_field(nu_restart,0,swvdf,'ruf8', &
           'swvdf',1,diag,field_loc_center, field_type_scalar)
      call read_restart_field(nu_restart,0,swidr,'ruf8', &
           'swidr',1,diag,field_loc_center, field_type_scalar)
      call read_restart_field(nu_restart,0,swidf,'ruf8', &
           'swidf',1,diag,field_loc_center, field_type_scalar)

      !-----------------------------------------------------------------
      ! ocean stress
      !-----------------------------------------------------------------
      if (my_task == master_task) &
           write(nu_diag,*) 'min/max ocean stress components'

      call read_restart_field(nu_restart,0,strocnxT,'ruf8', &
           'strocnxT',1,diag,field_loc_center, field_type_vector)
      call read_restart_field(nu_restart,0,strocnyT,'ruf8', &
           'strocnyT',1,diag,field_loc_center, field_type_vector)

      !-----------------------------------------------------------------
      ! internal stress
      ! The stress tensor must be read and scattered in pairs in order
      ! to properly match corner values across a tripole grid cut.
      !-----------------------------------------------------------------
      if (my_task == master_task) write(nu_diag,*) &
           'internal stress components'
      
      call read_restart_field(nu_restart,0,stressp_1,'ruf8', &
           'stressp_1',1,diag,field_loc_center,field_type_scalar) ! stressp_1
      call read_restart_field(nu_restart,0,stressp_3,'ruf8', &
           'stressp_3',1,diag,field_loc_center,field_type_scalar) ! stressp_3
      call read_restart_field(nu_restart,0,stressp_2,'ruf8', &
           'stressp_2',1,diag,field_loc_center,field_type_scalar) ! stressp_2
      call read_restart_field(nu_restart,0,stressp_4,'ruf8', &
           'stressp_4',1,diag,field_loc_center,field_type_scalar) ! stressp_4

      call read_restart_field(nu_restart,0,stressm_1,'ruf8', &
           'stressm_1',1,diag,field_loc_center,field_type_scalar) ! stressm_1
      call read_restart_field(nu_restart,0,stressm_3,'ruf8', &
           'stressm_3',1,diag,field_loc_center,field_type_scalar) ! stressm_3
      call read_restart_field(nu_restart,0,stressm_2,'ruf8', &
           'stressm_2',1,diag,field_loc_center,field_type_scalar) ! stressm_2
      call read_restart_field(nu_restart,0,stressm_4,'ruf8', &
           'stressm_4',1,diag,field_loc_center,field_type_scalar) ! stressm_4

      call read_restart_field(nu_restart,0,stress12_1,'ruf8', &
           'stress12_1',1,diag,field_loc_center,field_type_scalar) ! stress12_1
      call read_restart_field(nu_restart,0,stress12_3,'ruf8', &
           'stress12_3',1,diag,field_loc_center,field_type_scalar) ! stress12_1

      call read_restart_field(nu_restart,0,stress12_2,'ruf8', &
           'stress12_2',1,diag,field_loc_center,field_type_scalar) ! stress12_2
      call read_restart_field(nu_restart,0,stress12_4,'ruf8', &
           'stress12_4',1,diag,field_loc_center,field_type_scalar) ! stress12_4

      if (trim(grid_type) == 'tripole') then
         call ice_HaloUpdate_stress(stressp_1, stressp_3, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stressp_3, stressp_1, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stressp_2, stressp_4, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stressp_4, stressp_2, halo_info, &
                                    field_loc_center,  field_type_scalar)

         call ice_HaloUpdate_stress(stressm_1, stressm_3, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stressm_3, stressm_1, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stressm_2, stressm_4, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stressm_4, stressm_2, halo_info, &
                                    field_loc_center,  field_type_scalar)

         call ice_HaloUpdate_stress(stress12_1, stress12_3, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stress12_3, stress12_1, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stress12_2, stress12_4, halo_info, &
                                    field_loc_center,  field_type_scalar)
         call ice_HaloUpdate_stress(stress12_4, stress12_2, halo_info, &
                                    field_loc_center,  field_type_scalar)
      endif

      !-----------------------------------------------------------------
      ! ice mask for dynamics
      !-----------------------------------------------------------------
      if (my_task == master_task) &
           write(nu_diag,*) 'ice mask for dynamics'

      call read_restart_field(nu_restart,0,work1,'ruf8', &
           'iceumask',1,diag,field_loc_center, field_type_scalar)

      iceumask(:,:,:) = .false.
      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
      do iblk = 1, nblocks
         do j = 1, ny_block
         do i = 1, nx_block
            if (work1(i,j,iblk) > p5) iceumask(i,j,iblk) = .true.
         enddo
         enddo
      enddo
      !$OMP END PARALLEL DO

      ! for mixed layer model
      if (oceanmixed_ice) then

         if (my_task == master_task) &
              write(nu_diag,*) 'min/max sst, frzmlt'

         call read_restart_field(nu_restart,0,sst,'ruf8', &
              'sst',1,diag,field_loc_center, field_type_scalar)
         call read_restart_field(nu_restart,0,frzmlt,'ruf8', &
              'frzmlt',1,diag,field_loc_center, field_type_scalar)
      endif

      !-----------------------------------------------------------------
      ! Ensure unused stress values in west and south ghost cells are 0
      !-----------------------------------------------------------------
      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
      do iblk = 1, nblocks
         do j = 1, nghost
         do i = 1, nx_block
            stressp_1 (i,j,iblk) = c0
            stressp_2 (i,j,iblk) = c0
            stressp_3 (i,j,iblk) = c0
            stressp_4 (i,j,iblk) = c0
            stressm_1 (i,j,iblk) = c0
            stressm_2 (i,j,iblk) = c0
            stressm_3 (i,j,iblk) = c0
            stressm_4 (i,j,iblk) = c0
            stress12_1(i,j,iblk) = c0
            stress12_2(i,j,iblk) = c0
            stress12_3(i,j,iblk) = c0
            stress12_4(i,j,iblk) = c0
         enddo
         enddo
         do j = 1, ny_block
         do i = 1, nghost
            stressp_1 (i,j,iblk) = c0
            stressp_2 (i,j,iblk) = c0
            stressp_3 (i,j,iblk) = c0
            stressp_4 (i,j,iblk) = c0
            stressm_1 (i,j,iblk) = c0
            stressm_2 (i,j,iblk) = c0
            stressm_3 (i,j,iblk) = c0
            stressm_4 (i,j,iblk) = c0
            stress12_1(i,j,iblk) = c0
            stress12_2(i,j,iblk) = c0
            stress12_3(i,j,iblk) = c0
            stress12_4(i,j,iblk) = c0
         enddo
         enddo
      enddo
      !$OMP END PARALLEL DO

      !-----------------------------------------------------------------
      ! Ensure ice is binned in correct categories
      ! (should not be necessary unless restarting from a run with
      !  different category boundaries).
      !
      ! If called, this subroutine does not give exact restart.
      !-----------------------------------------------------------------
!!!      call cleanup_itd

      !-----------------------------------------------------------------
      ! compute aggregate ice state and open water area
      !-----------------------------------------------------------------

      !$OMP PARALLEL DO PRIVATE(iblk)
      do iblk = 1, nblocks

         call aggregate (nx_block, ny_block, &
                         aicen(:,:,:,iblk),  &
                         trcrn(:,:,:,:,iblk),&
                         vicen(:,:,:,iblk),  &
                         vsnon(:,:,:,iblk),  &
                         aice (:,:,  iblk),  &
                         trcr (:,:,:,iblk),  &
                         vice (:,:,  iblk),  &
                         vsno (:,:,  iblk),  &
                         aice0(:,:,  iblk),  &
                         tmask(:,:,  iblk),  &
                         max_ntrcr,          &
                         trcr_depend)

         aice_init(:,:,iblk) = aice(:,:,iblk)

      enddo
      !$OMP END PARALLEL DO

      ! if runid is bering then need to correct npt for istep0
      if (trim(runid) == 'bering') then
         npt = npt - istep0
      endif

      end subroutine restartfile

!=======================================================================

! Restarts from a CICE v4.1 (binary) dump
! author Elizabeth C. Hunke, LANL

      subroutine restartfile_v4 (ice_ic)

      use ice_broadcast, only: broadcast_scalar
      use ice_blocks, only: nghost, nx_block, ny_block
      use ice_calendar, only: istep0, istep1, time, time_forc, calendar, npt
      use ice_communicate, only: my_task, master_task
      use ice_constants, only: c0, p5, &
          field_loc_center, field_loc_NEcorner, &
          field_type_scalar, field_type_vector
      use ice_domain, only: nblocks, distrb_info
      use ice_domain_size, only: nilyr, nslyr, ncat, nx_global, ny_global, &
          max_ntrcr, max_blocks
      use ice_fileunits, only: nu_diag, nu_rst_pointer, nu_restart
      use ice_flux, only: scale_factor, swvdr, swvdf, swidr, swidf, &
          strocnxT, strocnyT, sst, frzmlt, iceumask, &
          stressp_1, stressp_2, stressp_3, stressp_4, &
          stressm_1, stressm_2, stressm_3, stressm_4, &
          stress12_1, stress12_2, stress12_3, stress12_4
      use ice_gather_scatter, only: scatter_global_stress
      use ice_grid, only: tmask
      use ice_itd, only: aggregate
      use ice_ocean, only: oceanmixed_ice
      use ice_read_write, only: ice_open, ice_read, ice_read_global
      use ice_state, only: trcr_depend, aice, vice, vsno, trcr, &
          aice0, aicen, vicen, vsnon, trcrn, aice_init, &
          nt_Tsfc, nt_sice, nt_qice, nt_qsno, uvel, vvel

      character (*), optional :: ice_ic

      ! local variables

      integer (kind=int_kind) :: &
         i, j, k, n, iblk, &     ! counting indices
         iignore                 ! dummy variable

      real (kind=real_kind) :: &
         rignore                 ! dummy variable

      character(len=char_len_long) :: &
         filename, filename0

      logical (kind=log_kind) :: &
         diag

      real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
         work1

      real (kind=dbl_kind), dimension(:,:), allocatable :: &
         work_g1, work_g2

      if (present(ice_ic)) then
         filename = ice_ic
      elseif (my_task == master_task) then
         open(nu_rst_pointer,file=pointer_file)
         read(nu_rst_pointer,'(a)') filename0
         filename = trim(filename0)
         close(nu_rst_pointer)
         write(nu_diag,*) 'Read ',pointer_file(1:lenstr(pointer_file))
      endif

      call ice_open(nu_restart,filename,0)

      if (my_task == master_task) &
         write(nu_diag,*) 'Using restart dump=', trim(filename)

      if (use_restart_time) then

         if (my_task == master_task) then
            read (nu_restart) istep0,time,time_forc
            write(nu_diag,*) 'Restart read at istep=',istep0,time,time_forc
         endif
         call broadcast_scalar(istep0,master_task)
         istep1 = istep0
         call broadcast_scalar(time,master_task)
         call broadcast_scalar(time_forc,master_task)
         call calendar(time)

      else

         if (my_task == master_task) &
            read (nu_restart) iignore,rignore,rignore

      endif

      diag = .true.     ! write min/max diagnostics for field

      !-----------------------------------------------------------------
      ! state variables
      ! Tsfc is the only tracer read in this file.  All other
      ! tracers are in their own dump/restart files.
      !-----------------------------------------------------------------
      do n=1,ncat
         if (my_task == master_task) &
              write(nu_diag,*) 'cat ',n, &
                               ' min/max area, vol ice, vol snow, Tsfc'

         call ice_read(nu_restart,0,aicen(:,:,n,:),'ruf8',diag, &
                       field_loc_center, field_type_scalar)
         call ice_read(nu_restart,0,vicen(:,:,n,:),'ruf8',diag, &
                       field_loc_center, field_type_scalar)
         call ice_read(nu_restart,0,vsnon(:,:,n,:),'ruf8',diag, &
                       field_loc_center, field_type_scalar)
         call ice_read(nu_restart,0,trcrn(:,:,nt_Tsfc,n,:),'ruf8',diag, &
                       field_loc_center, field_type_scalar)

         if (my_task == master_task) &
              write(nu_diag,*) 'cat ',n, 'min/max sice for each layer'
         do k=1,nilyr
            call ice_read(nu_restart,0,trcrn(:,:,nt_sice+k-1,n,:),'ruf8',diag, &
                          field_loc_center, field_type_scalar)
         enddo

         if (my_task == master_task) &
              write(nu_diag,*) 'cat ',n, 'min/max qice for each layer'
         do k=1,nilyr
            call ice_read(nu_restart,0,trcrn(:,:,nt_qice+k-1,n,:),'ruf8',diag, &
                          field_loc_center, field_type_scalar)
         enddo

         if (my_task == master_task) &
              write(nu_diag,*) 'cat ',n, 'min/max qsno for each layer'
         do k=1,nslyr
            call ice_read(nu_restart,0,trcrn(:,:,nt_qsno+k-1,n,:),'ruf8',diag, &
                          field_loc_center, field_type_scalar)
         enddo
      enddo ! ncat

      !-----------------------------------------------------------------
      ! velocity
      !-----------------------------------------------------------------
      if (my_task == master_task) &
           write(nu_diag,*) 'min/max velocity components'

      call ice_read(nu_restart,0,uvel,'ruf8',diag, &
                       field_loc_NEcorner, field_type_vector)
      call ice_read(nu_restart,0,vvel,'ruf8',diag, &
                       field_loc_NEcorner, field_type_vector)

      !-----------------------------------------------------------------
      ! radiation fields
      !-----------------------------------------------------------------

      if (my_task == master_task) &
         write(nu_diag,*) 'radiation fields'

      call ice_read(nu_restart,0,scale_factor,'ruf8',diag, &
                    field_loc_center, field_type_scalar)
      call ice_read(nu_restart,0,swvdr,'ruf8',diag, &
                    field_loc_center, field_type_scalar)
      call ice_read(nu_restart,0,swvdf,'ruf8',diag, &
                    field_loc_center, field_type_scalar)
      call ice_read(nu_restart,0,swidr,'ruf8',diag, &
                    field_loc_center, field_type_scalar)
      call ice_read(nu_restart,0,swidf,'ruf8',diag, &
                    field_loc_center, field_type_scalar)

      !-----------------------------------------------------------------
      ! ocean stress
      !-----------------------------------------------------------------
      if (my_task == master_task) &
           write(nu_diag,*) 'min/max ocean stress components'

      call ice_read(nu_restart,0,strocnxT,'ruf8',diag, &
                       field_loc_center, field_type_vector)
      call ice_read(nu_restart,0,strocnyT,'ruf8',diag, &
                       field_loc_center, field_type_vector)

      !-----------------------------------------------------------------
      ! internal stress
      ! The stress tensor must be read and scattered in pairs in order
      ! to properly match corner values across a tripole grid cut.
      !-----------------------------------------------------------------
      if (my_task == master_task) write(nu_diag,*) &
           'internal stress components'
      
      allocate (work_g1(nx_global,ny_global), &
                work_g2(nx_global,ny_global))

      call ice_read_global(nu_restart,0,work_g1,'ruf8',diag) ! stressp_1
      call ice_read_global(nu_restart,0,work_g2,'ruf8',diag) ! stressp_3
      call scatter_global_stress(stressp_1, work_g1, work_g2, &
                                 master_task, distrb_info)
      call scatter_global_stress(stressp_3, work_g2, work_g1, &
                                 master_task, distrb_info)

      call ice_read_global(nu_restart,0,work_g1,'ruf8',diag) ! stressp_2
      call ice_read_global(nu_restart,0,work_g2,'ruf8',diag) ! stressp_4
      call scatter_global_stress(stressp_2, work_g1, work_g2, &
                                 master_task, distrb_info)
      call scatter_global_stress(stressp_4, work_g2, work_g1, &
                                 master_task, distrb_info)

      call ice_read_global(nu_restart,0,work_g1,'ruf8',diag) ! stressm_1
      call ice_read_global(nu_restart,0,work_g2,'ruf8',diag) ! stressm_3
      call scatter_global_stress(stressm_1, work_g1, work_g2, &
                                 master_task, distrb_info)
      call scatter_global_stress(stressm_3, work_g2, work_g1, &
                                 master_task, distrb_info)

      call ice_read_global(nu_restart,0,work_g1,'ruf8',diag) ! stressm_2
      call ice_read_global(nu_restart,0,work_g2,'ruf8',diag) ! stressm_4
      call scatter_global_stress(stressm_2, work_g1, work_g2, &
                                 master_task, distrb_info)
      call scatter_global_stress(stressm_4, work_g2, work_g1, &
                                 master_task, distrb_info)

      call ice_read_global(nu_restart,0,work_g1,'ruf8',diag) ! stress12_1
      call ice_read_global(nu_restart,0,work_g2,'ruf8',diag) ! stress12_3
      call scatter_global_stress(stress12_1, work_g1, work_g2, &
                                 master_task, distrb_info)
      call scatter_global_stress(stress12_3, work_g2, work_g1, &
                                 master_task, distrb_info)

      call ice_read_global(nu_restart,0,work_g1,'ruf8',diag) ! stress12_2
      call ice_read_global(nu_restart,0,work_g2,'ruf8',diag) ! stress12_4
      call scatter_global_stress(stress12_2, work_g1, work_g2, &
                                 master_task, distrb_info)
      call scatter_global_stress(stress12_4, work_g2, work_g1, &
                                 master_task, distrb_info)

      deallocate (work_g1, work_g2)

      !-----------------------------------------------------------------
      ! ice mask for dynamics
      !-----------------------------------------------------------------
      if (my_task == master_task) &
           write(nu_diag,*) 'ice mask for dynamics'

      call ice_read(nu_restart,0,work1,'ruf8',diag, &
                    field_loc_center, field_type_scalar)

      iceumask(:,:,:) = .false.
      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
      do iblk = 1, nblocks
         do j = 1, ny_block
         do i = 1, nx_block
            if (work1(i,j,iblk) > p5) iceumask(i,j,iblk) = .true.
         enddo
         enddo
      enddo
      !$OMP END PARALLEL DO

      ! for mixed layer model
      if (oceanmixed_ice) then

         if (my_task == master_task) &
              write(nu_diag,*) 'min/max sst, frzmlt'

         call ice_read(nu_restart,0,sst,'ruf8',diag, &
                       field_loc_center, field_type_scalar)
         call ice_read(nu_restart,0,frzmlt,'ruf8',diag, &
                       field_loc_center, field_type_scalar)
      endif

      if (my_task == master_task) close(nu_restart)

      !-----------------------------------------------------------------
      ! Ensure unused stress values in west and south ghost cells are 0
      !-----------------------------------------------------------------
      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
      do iblk = 1, nblocks
         do j = 1, nghost
         do i = 1, nx_block
            stressp_1 (i,j,iblk) = c0
            stressp_2 (i,j,iblk) = c0
            stressp_3 (i,j,iblk) = c0
            stressp_4 (i,j,iblk) = c0
            stressm_1 (i,j,iblk) = c0
            stressm_2 (i,j,iblk) = c0
            stressm_3 (i,j,iblk) = c0
            stressm_4 (i,j,iblk) = c0
            stress12_1(i,j,iblk) = c0
            stress12_2(i,j,iblk) = c0
            stress12_3(i,j,iblk) = c0
            stress12_4(i,j,iblk) = c0
         enddo
         enddo
         do j = 1, ny_block
         do i = 1, nghost
            stressp_1 (i,j,iblk) = c0
            stressp_2 (i,j,iblk) = c0
            stressp_3 (i,j,iblk) = c0
            stressp_4 (i,j,iblk) = c0
            stressm_1 (i,j,iblk) = c0
            stressm_2 (i,j,iblk) = c0
            stressm_3 (i,j,iblk) = c0
            stressm_4 (i,j,iblk) = c0
            stress12_1(i,j,iblk) = c0
            stress12_2(i,j,iblk) = c0
            stress12_3(i,j,iblk) = c0
            stress12_4(i,j,iblk) = c0
         enddo
         enddo
      enddo
      !$OMP END PARALLEL DO

      !-----------------------------------------------------------------
      ! Ensure ice is binned in correct categories
      ! (should not be necessary unless restarting from a run with
      !  different category boundaries).
      !
      ! If called, this subroutine does not give exact restart.
      !-----------------------------------------------------------------
!!!      call cleanup_itd

      !-----------------------------------------------------------------
      ! compute aggregate ice state and open water area
      !-----------------------------------------------------------------

      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
      do iblk = 1, nblocks

         call aggregate (nx_block, ny_block, &
                         aicen(:,:,:,iblk),  &
                         trcrn(:,:,:,:,iblk),&
                         vicen(:,:,:,iblk),  &
                         vsnon(:,:,:,iblk),  &
                         aice (:,:,  iblk),  &
                         trcr (:,:,:,iblk),  &
                         vice (:,:,  iblk),  &
                         vsno (:,:,  iblk),  &
                         aice0(:,:,  iblk),  &
                         tmask(:,:,  iblk),  &
                         max_ntrcr,          &
                         trcr_depend)

         aice_init(:,:,iblk) = aice(:,:,iblk)

      enddo
      !$OMP END PARALLEL DO

      ! creates netcdf if restart_format = 'nc'
      filename = trim(restart_dir) // '/iced.converted'
      call dumpfile(filename) 
      call final_restart
      ! stop

      ! if runid is bering then need to correct npt for istep0
      if (trim(runid) == 'bering') then
         npt = npt - istep0
      endif

      end subroutine restartfile_v4

!=======================================================================

      end module ice_restart_driver

!=======================================================================
