module m_evp_restart

#if defined (EVP_IO_TAUICE)
      integer, private,parameter :: nfields=4
#else
      integer, private,parameter :: nfields=2
#endif

      contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


! Routine to read EVP fields

      integer function evp_read_restart(rt,nstep0,dtime0)
      use mod_evp
      use mod_year_info
      use mod_za
      use mod_common_ice, only: tauxice, tauyice
      use mod_hycom_nersc
      implicit none
      type(year_info)     , intent(in)  :: rt
      integer,              intent(in)  :: nstep0
      real*8,               intent(in)  :: dtime0
      character(len=50) :: filename
      integer :: j,iostat
      logical :: ex



      character*80    l_ctitle(4)
      integer :: l_iversn, l_iexpt, l_yrflag,l_itdm, 
     &           l_jtdm, iyy, idd, ihh, iss

      character(len=50) :: filebase
      character(len=8) :: char8
      real :: xmin, xmax
      integer :: nop,ifld

      real :: xmin2, xmax2, l_time
      integer :: l_nstep,l_imem,offset,k,offsetb,ios,ifldb,nfields_input
      logical :: first

#if defined (SINGLE_RESTART) && defined (SINGLE_RESTART_ONLY)
      ! In this case this routine is entirely replaced by the following
      ! call
      evp_read_restart=evp_read_restart_new(rt)
      return
#endif

      ! New output format -- uses hycom-style IO (zaio)
      nop=710
      filebase=rungen//'restart'//rt%cyy//'_'//rt%cdd//'_'//rt%chh//
     &         'EVP'
      inquire (exist=ex, file=trim(filebase)//'.a')
      if (.not.ex) then
         if (mnproc==1) write(lp,*)'No such file '//trim(filebase)//'.a'
         evp_read_restart=1
         return
      else
         if (mnproc==1) then
            write(lp,'(a,i4)')'evp_read_restart: reading '
     &         //trim(filebase)//'.[ab] record=',imem
         end if
      end if

      call zaiopf(trim(filebase)//'.a','old',nop)
      open(nop,file=trim(filebase)//'.b',status='old',
     &         form='formatted',access='direct',recl=80)

      ! Number of fields for one member (only u/v)
      !KAL - TODO - must fix the offset/nfields logic when introducing new fields

      ! Offset for IO
      ! KAL -- offset not used anymore on input - search for correct
      ! KAL -- imem on input in stead
      offset=(imem-1)*nfields
      offsetb=offset+1 ! First record in .b  is model info

      ! Read "header"
      read(nop,118,rec=1) l_itdm,l_jtdm
      if (itdm/=l_itdm .or. jtdm /= l_jtdm) then
         if (mnproc==1) write(lp,*) 'Dimension mismatch in EVP file'
         call xcstop('(evp_read_restart)')
         stop '(evp_read_restart)'
      end if


      ! Read .b file until we find correct imem
      first=.true.
      ios=0
      ifldb=1 ! First record was read above
      nfields_input=0
      do while (ios==0) 

         read(nop,117,rec=ifldb+1,iostat=ios)
     &      char8,l_nstep,l_time,l_imem,xmin2,xmax2
         ifldb=ifldb+1

         ! Locate first field belonging to imem
         if (first) then
            if (l_imem==imem .and. ios==0) then
               offsetb=ifldb
               first=.false.
            end if
         end if

         ! Update num fields belonging to imem
         if (l_imem==imem .and. ios==0) nfields_input=nfields_input+1
         !print *,ifldb,nfields_input,ios
      end do

      ! Here we never found member imem
      if (first) then
         if (mnproc==1) write(lp,'(a,i3)') 'Could not find member ',imem
         evp_read_restart=1
         return
      end if

      ! Offset into ".a" - file
      offset=offsetb-1

      ! Adjust offset so that offset+1 is first record
      offset =offset -1
      offsetb=offsetb-1

      ! Info 
      if (mnproc==1) then
         write(lp,'(a,i4,a,i4)')   'evp_read_restart: restart for ',
     &      imem,' starts at .a offset ',offset+1
         write(lp,'(a,i4,a,i4,a)') 'evp_read_restart: restart for ',
     &      imem,' contains ',nfields_input,' model fields'
      end if
            


      ! Skip offset number of records in .a-file
      do k=1,offset
         call zaiosk(nop)
      end do


      ! Start reading at correct offsets
      do ifld=1,nfields_input

         ! Same order as out...
         read(nop,117,rec=offsetb+ifld,iostat=ios)
     &      char8,l_nstep,l_time,l_imem,xmin2,xmax2

         if (l_imem/=imem) then
            if (mnproc==1) write(lp,*) 'Ensemble member mismatch..'
            call xcstop('(evp_read_restart)')
            stop '(evp_read_restart)'

            ! TODO Change this with "cycle" ?
         end if

         if (mnproc==1)
     &   write(lp,*)char8,imem,xmin2,xmax2

         if     (trim(char8)=='uvel') then
            call zaiord(uvel(1-nbdy,1-nbdy),ip,.false.,
     &                  xmin,xmax, nop)
         elseif (trim(char8)=='vvel') then
            call zaiord(vvel(1-nbdy,1-nbdy),ip,.false.,
     &                  xmin,xmax, nop)
         elseif     (trim(char8)=='tauxice') then
            call zaiord(tauxice(1-nbdy,1-nbdy),ip,.false.,
     &                  xmin,xmax, nop)
         elseif (trim(char8)=='tauyice') then
            call zaiord(tauyice(1-nbdy,1-nbdy),ip,.false.,
     &                  xmin,xmax, nop)
         else
            if (mnproc==1) print *, 'Unknown field '//char8
            call xcstop('(evp_read_restart')
            stop '(evp_read_restart'
         end if

         ! Check for coherency between min and maxes
         if (abs(xmin-xmin2)>abs(xmin)*1e-4 .or.
     &       abs(xmax-xmax2)>abs(xmax)*1e-4) then
            if     (mnproc.eq.1) then
               write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &           'error - .a and .b files not consistent:',
     &           '.a,.b min = ',xmin2,xmin ,xmin2-xmin ,
     &           '.a,.b max = ',xmax2,xmax ,xmax2-xmax
               print *,'Variable :',char8
            endif
            call xcstop('(evp_read_restart)')
            stop '(evp_read_restart)'
         end if


      end do

      call zaiocl(nop)
      close(nop)
      evp_read_restart=0


 117  format (a8,' =',i11,f11.2,i3,1p2e16.7)
 118  format ('itdm= ',i5, ' jtdm= ',i5)
      end function evp_read_restart



! Routine to save EVP fields

      integer function evp_save_restart(rt,nstep0,dtime0)
      use mod_xc
      use mod_za
      use mod_evp
      use mod_year_info
      use mod_hycom_nersc
      use mod_common_ice, only: tauxice, tauyice
      implicit none
      type(year_info)     , intent(in) :: rt
      integer,              intent(in) :: nstep0
      real*8,               intent(in) :: dtime0
      character(len=50) :: filename
      integer :: j,iostat,ios
      integer ioproc, ierr
      
      integer :: nop,offset,k,l_nstep,l_imem,offsetb
      real :: xmin,xmax, l_time
      character(len=50) :: filebase
      character(len= 1) :: char1
      character(len= 8) :: char8


#if defined (SINGLE_RESTART) /* Dump to new restart system */
      evp_save_restart=evp_save_restart_new(rt)
#if defined (SINGLE_RESTART_ONLY) /* This short circuits this return */
      return
#endif
#endif

      evp_save_restart=0

      ! New output format -- uses hycom-style IO (zaio)
      nop=710
      filebase=rungen//'restart'//rt%cyy//'_'//
     &         rt%cdd//'_'//rt%chh//'EVP'
      if (mnproc==1) then
         write(lp,'(a,i4)',advance='yes') 'evp_save_restart: saving '//
     &         trim(filebase)//'.[ab] record=',imem
      end if
      call zaiopf(trim(filebase)//'.a','unknown',nop)
      if (mnproc==1) then
         open(nop,file=trim(filebase)//'.b',status='unknown',
     &        form='formatted',access='direct',recl=80)
      end if


      ! Offset for IO
      offset=(imem-1)*nfields
      offsetb=offset+1 ! First record is model info 
      if (mnproc==1) then
         write(nop,118,rec=1) itdm,jtdm
         call flush(nop)
      end if

      ! Skip offset number of records
      do k=1,offset
         call zaiosk(nop)
      end do


      ! The rest ...
      call zaiowr(uvel(1-nbdy,1-nbdy),ip,.false.,
     &            xmin,xmax, nop, .false.)
      if (mnproc==1) then
         write(nop,117,rec=offsetb+1)
     &      'uvel    ',nstep0,dtime0,imem,xmin,xmax
         call flush(nop)
      end if

      call zaiowr(vvel(1-nbdy,1-nbdy),ip,.false.,
     &            xmin,xmax, nop, .false.)
      if (mnproc==1) then
         write(nop,117,rec=offsetb+2)
     &      'vvel    ',nstep0,dtime0,imem,xmin,xmax
         call flush(nop)
      end if


#if defined (EVP_IO_TAUICE)
      call zaiowr(tauxice(1-nbdy,1-nbdy),ip,.false.,
     &            xmin,xmax, nop, .false.)
      if (mnproc==1) then
         write(nop,117,rec=offsetb+3)
     &      'tauxice ',nstep0,dtime0,imem,xmin,xmax
         call flush(nop)
      end if

      call zaiowr(tauyice(1-nbdy,1-nbdy),ip,.false.,
     &            xmin,xmax, nop, .false.)
      if (mnproc==1) then
         write(nop,117,rec=offsetb+4)
     &      'tauyice ',nstep0,dtime0,imem,xmin,xmax
         call flush(nop)
      end if
#endif

      call zaiocl(nop)
      if (mnproc==1)  close(nop)
      evp_save_restart=0

      return

 117  format (a8,' =',i11,f11.2,i3,1p2e16.7)
 118  format ('itdm= ',i5, ' jtdm= ',i5)


      end function evp_save_restart




#if defined (SINGLE_RESTART)
      integer function evp_save_restart_new(rt)
      use mod_evp
      use mod_year_info
      use mod_za
      use mod_common_ice, only: tauxice, tauyice
      use mod_hycom_nersc
      use mod_restart
      implicit none
      type(year_info)     , intent(in) :: rt
      
      integer :: ierr(4)
      character(len=80) :: filebase

      ierr=0
      evp_save_restart_new=0
      filebase=trim(restartfile(rt))
      if (mnproc==1) then
         write(lp,'(a)')'evp_save_restart_new: reading '
     &      //trim(filebase)//'.[ab]'
      end if

      call restart_appendfield(filebase,'uice    ',0,1,uvel    ,ierr(1))
      call restart_appendfield(filebase,'vice    ',0,1,vvel    ,ierr(2))
      call restart_appendfield(filebase,'tauxice ',0,1,tauxice ,ierr(3))
      call restart_appendfield(filebase,'tauyice ',0,1,tauyice ,ierr(4))
      evp_save_restart_new=maxval(abs(ierr))
      end function evp_save_restart_new


      function evp_read_restart_new(rt)
      use mod_xc
      use mod_za
      use mod_evp
      use mod_common_ice, only : tauxice, tauyice
      use mod_year_info
      use mod_hycom_nersc
      use mod_restart
      implicit none
      type(year_info) , intent(in) :: rt

      character filebase*80
      integer :: evp_read_restart_new
      integer :: ierr(4)

      evp_read_restart_new=0
      filebase=trim(restartfile(rt))
      if (mnproc==1) then
         write(lp,'(a)')'evp_read_restart_new: reading '
     &      //trim(filebase)//'.[ab]'
      end if

      ! Retrieve record indexes for the variables
      call restart_getfield(filebase,'uice'   ,0,1,uvel   ,ierr(1))
      call restart_getfield(filebase,'vice'   ,0,1,vvel   ,ierr(2))
      call restart_getfield(filebase,'tauxice',0,1,tauxice,ierr(3))
      call restart_getfield(filebase,'tauyice',0,1,tauyice,ierr(4))

      ! These are not essential fields so we warn if they are not found
      if (any(ierr(1:4)/=0)) then
         if (mnproc==1) write (lp,'(a)') 'An error occured when '// 
     &      ' reading the evp fields'
         evp_read_restart_new=-1
      end if
      end function evp_read_restart_new
#endif
      end module m_evp_restart