MODULE m_read_icemod
! Francois Massonnet, UCL, 2013
! Reads data from icemod file (instead of classically restart files).
! This is required when doing data assimilation of OSISAF sea ice drift
! computed over several hours/days. Indeed, the restart only gives a
! snapshot of the state of the system while the icemod records the time
! average. The icemod file should have one time slice.
USE NETCDF
#if defined (QMPI)
  use qmpi
#else
  use qmpi_fake
#endif

CONTAINS
        SUBROUTINE read_icemod(fld,k,enslist,cfld,nx,ny)
        IMPLICIT NONE

        real,dimension(nx,ny),intent(inout):: fld    ! output fl
        character(len=*),       intent(in) :: cfld   ! name of fld
        integer,                intent(in) :: k      ! Index to enslist
        integer,dimension(:),   intent(in) :: enslist! List of existing ensemble members
        integer,                intent(in) :: nx,ny  ! Grid dimension
        integer                                 :: iens
        integer                         :: error, ncid,varID
        character(len=3)  :: cmem
        character(len=99) :: cfile
        logical           :: exf
        iens = enslist(k)
        write(cmem,'(i3.3)') 100+iens  ! iens=1 gives cmem = 101
        
        cfile='icemod_'//cmem//'.nc'

        inquire(file=cfile,  exist=exf)
        if (.not.exf) then
                if (master) print *, '(read_icemod): Icemod file '//cfile//' missing!'
                call stop_mpi()
        end if

        error = nf90_open(trim(cfile),nf90_Write,ncid); if (error.ne.nf90_noerr) call handle_err(error, "opening")
        error = nf90_inq_varid(ncid, trim(cfld), varID); if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")       
        error = nf90_get_var(ncid, varID, fld); if (error.ne.nf90_noerr) call handle_err(error, "getting 2D variable")
 
        END SUBROUTINE read_icemod

        subroutine handle_err(status, infomsg)
          integer,            intent ( in) :: status
          character(len = *), intent ( in), optional :: infomsg
          if(status /= nf90_noerr) then
             if (master) then
                if (present(infomsg)) then
                   print *, 'Error while '//infomsg//' - '//trim(nf90_strerror(status))
                else
                   print *, trim(nf90_strerror(status))
                endif ! opt arg
                print *,'(io_mod_fld)'
             endif ! only master outputs
             call stop_mpi()
          end if ! check error status
        end subroutine handle_err
 

END MODULE m_read_icemod