

c
c   program to convert restart files
c   Here not to netcdf, but to include e.g. more tracers
c   or to convert from single to double precision
c

c#define change_time_stamp
c#define  add_ledwell_tracer
c#define with_ktmix
#define with_tkemix

c#define convert_real4_to_real8
#define convert_real8_to_real4

      program cv_rest
      implicit none
      integer :: io,nt,itt,ierr, io2,nt2,io3
      character (len=32) :: stamp
      character (len=40) :: tr_name

      integer :: imt,jmt,km  , i,j,k,n
      integer, allocatable :: kmt(:,:)

#ifdef convert_real4_to_real8
      real (kind=4), allocatable  :: dat_in(:,:)
      real (kind=8), allocatable  :: dat_out(:,:)
#elif defined convert_real8_to_real4
      real (kind=8), allocatable  :: dat_in(:,:)
      real (kind=4), allocatable  :: dat_out(:,:)
#else
      real, allocatable  :: dat_in(:,:)
      real, allocatable  :: dat_out(:,:)
#endif

      real, allocatable  :: tr(:,:,:,:)
      real :: x
      integer :: indp
      real,allocatable :: xt(:), yt(:), xu(:), yu(:), zt(:),zw(:)

      print*,' Converting '

      call get_commandline_arg_int(1,imt,ierr)
      if (ierr/=0) goto 99
      call get_commandline_arg_int(2,jmt,ierr)
      if (ierr/=0) goto 99
      call get_commandline_arg_int(3,km,ierr)
      if (ierr/=0) goto 99
      call get_commandline_arg_int(4,nt2,ierr)
      if (ierr/=0) goto 99

      allocate(kmt(imt,jmt) )
      allocate(dat_in(imt,jmt), dat_out(imt,jmt) )
c
c     read topography to zero out special values
c
      call getunit (io, 'kmt.dta', 'usr ieee')
      read (io) ! no checks
c       skip header record and read dimensional info
      read (io) 
      read (io) i, j, k
      read (io) 
      read (io) kmt
      close(io)

#ifdef add_ledwell_tracer

      print*,
     & ' =>Reading the grid definition from file grid.dta'
      call getunit (io, 'grid.dta', 'u s r ieee')
      read (io) ! no checks for awi id
      read (io) i, j, k
      allocate( xt(imt), yt(jmt), xu(imt), yu(jmt))
      allocate( zt(km),zw(km) )
      read (io) (x,i=1,imt)
     &,         (x,j=1,jmt)
     &,         (x,i=1,imt)
     &,         (x,j=1,jmt)
     &,         (x,k=1,km)
     &,         (x,k=0,km)
     &,         (xt(i),i=1,imt)
     &,         (xu(i),i=1,imt)
     &,         (yt(j),j=1,jmt)
     &,         (yu(j),j=1,jmt)
     &,         (zt(k),k=1,km)
     &,         (zw(k),k=1,km)
      close(io)
#endif


      call get_free_iounit(io)
      open(io,file='old_restart.dta',form='unformatted',status='old')

      read(io) itt,stamp,nt

      call get_free_iounit(io2)
      open(io2,file='new_restart.dta',form='unformatted',status='new')

      print*,' time stamp = ',stamp

#ifdef change_time_stamp
c      stamp ='m/d/y= 2/ 1/1990, h:m:s=22:30: 0' 
      stamp ='m/d/y= 1/14/1990, h:m:s=16:30:00' 
      print*,' new time stamp = ',stamp
#endif


      print*,' itt        = ',itt
      print*,' nt         = ',nt
      print*,' nt2        = ',nt2
      write(io2) itt,stamp,nt2

      allocate(tr(imt,jmt,km,nt+1:nt2))

#ifdef add_ledwell_tracer
      tr=0.

c      i=indp(360.0-(80.2),xt,imt)
c      j=indp(28.0,yt,jmt)
c      k=indp(150.0e2,zt,km)
c      print*,' adding a dye at i,j,k=',i,j,k,' x,y,z=',xt(i),yt(j),zt(k)
c      print*,' for tracer # 3'
c      tr(i,j,k,3)=1.0

      i=indp(360.0-(28.0+15./60.),xt,imt)
      j=indp(25.0+40.0/60.0,yt,jmt)
      k=indp(310.0e2,zt,km)
      print*,' adding a dye at i,j,k=',i,j,k,' x,y,z=',xt(i),yt(j),zt(k)
      print*,' for tracer # 3'
      tr(i,j,k,3)=1.0

      i=indp(360.0-40.0,xt,imt)
      j=indp(20.0,yt,jmt)
      k=indp(332.7e2,zt,km)
      print*,' adding a dye at i,j,k=',i,j,k,' x,y,z=',xt(i),yt(j),zt(k)
      print*,' for tracer # 4'
      tr(i,j,k,4)=1.0

      i=indp(360.0-35.0,xt,imt)
      j=indp(22.0,yt,jmt)
      k=indp(332.7e2,zt,km)
      print*,' adding a dye at i,j,k=',i,j,k,' x,y,z=',xt(i),yt(j),zt(k)
      print*,' for tracer # 5'
      tr(i,j,k,5)=1.0

      i=indp(360.0-(60.),xt,imt)
      j=indp(42.5,yt,jmt)
      k=indp(1500.0e2,zt,km)
      print*,' adding a dye at i,j,k=',i,j,k,' x,y,z=',xt(i),yt(j),zt(k)
      print*,' for tracer # 3'
      tr(i,j,k,3)=1.0


      i=indp(360.0-76.0,xt,imt)
      j=indp(26.0,yt,jmt)
      k=indp(1438.7e2,zt,km)
      print*,' adding a dye at i,j,k=',i,j,k,' x,y,z=',xt(i),yt(j),zt(k)
      print*,' for tracer # 4'
      tr(i,j,k,4)=1.0

#else
c
c     read initial data
c
      do n=nt+1,nt2
        write(tr_name, '("tracer_",i2,".mom.ic")') n
        call replace_space_zero(tr_name)
        print*,' reading initial conditions from file ',
     &            tr_name(1:len_trim(tr_name))
        call getunit (io3,tr_name, 'u s r ieee')
        read(io3)  
        do j=1,jmt
         read(io3)  
         read(io3) stamp,x,k,k,k,x,k,(x,i=1,imt),(x,k=1,km),tr(:,j,:,n)
         do i=1,imt
          if (kmt(i,j)<km) tr(i,j,kmt(i,j)+1:km,n)=0.
         enddo
        enddo
        close(io3)
      enddo
#endif
c
c     barotropic mode
c
      print*,' reading writing barotropic mode '
      read(io)  dat_in ! ptd
      dat_out=dat_in
      write(io2)dat_out

      read(io)  dat_in ! guess
      dat_out=dat_in
      write(io2)dat_out

      read(io)  dat_in! psi(,,1)
      print*, 'psi(,,1) min,max [Sv]:',
     &                  minval(dat_in)/1.e12,maxval(dat_in)/1.e12
      dat_out=dat_in
      write(io2)dat_out

      read(io)  dat_in ! psi(,,2)
      print*, 'psi(,,2) min,max [Sv]:',
     &                  minval(dat_in)/1.e12,maxval(dat_in)/1.e12
      dat_out=dat_in
      write(io2)dat_out

      do k=1,km
       print*,' writing tracers k= ',k,' von ',km

       do n=1,nt ! read only nt_in tracers
        read(io) dat_in ! t(,k,,n,taum1)
        dat_out=dat_in
        write(io2) dat_out
        read(io) dat_in ! t(,k,,n,tau)
        dat_out=dat_in
        write(io2) dat_out
       enddo

c      add new tracer
       do n=nt+1,nt2
        write(io2) tr(:,:,k,n)
        write(io2) tr(:,:,k,n)
       enddo

       do n=1,2
        read(io) dat_in ! u(,k,,n,taum1)
        dat_out=dat_in
        write(io2) dat_out
        read(io) dat_in ! u(,k,,n,tau)
        dat_out=dat_in
        write(io2) dat_out
       enddo
      enddo

#ifdef with_ktmix
      ! ktmix ?
       print*,' writing dml '
      read(io) dat_in ! dml
      dat_out=dat_in
      write(io2) dat_out
#endif

#ifdef with_tkemix
      ! tkemix
      do k=1,km
       print*,' writing eke k= ',k,' von ',km
        read(io) dat_in ! eke(,k,,n,taum1)
        dat_out=dat_in
        write(io2) dat_out

        read(io) dat_in ! eke(,k,,n,tau)
        dat_out=dat_in
        write(io2) dat_out
      enddo
#endif

      close(io)
      close(io2)

      stop
 99   print*,' program needs command line parameter:'
      print*,'<executable> imt jmt km nt2'

      end program cv_rest



      subroutine halt_stop(s)
      character*(*) s
      print*,s
      stop
      end
