program redate_restart_mk3

! Reads a restart file written by the CSIRO Mk3 AGCM and displays the date.
! The user is prompted to enter a new date and a new restart file is
! written, which is identical to the first but with the new date.
!
! Usage:   redate_restart_mk3 <input_file> <output_file>
!
! Steve Phipps  5 December 2001

  implicit none

! Define parameters
  integer, parameter :: inunit=10, outunit=11, ngrid=13, ms=6

! Declare variables to hold data from restart file
  character(len=50) :: exptyp, header2
  integer :: ndays1, ndays2, mstep
  complex, dimension(:), allocatable :: spec
  real, dimension(:, :), allocatable :: qout, griddata, ron, son, stice, &
                                        ochfa, ssdnn, sgsav, oprec1
  integer, dimension(:, :), allocatable :: imsl, isflag
  real, dimension(:, :, :), allocatable :: sdot, wb, tggsn, qlb
  real, dimension(:, :, :, :), allocatable :: bigxt

! Declare other variables
  character(len=80) :: infile, outfile
  character(len=3) :: hres, vres
  integer :: i, n, nw, nl, lon, ln2, lat, lat2, lonc, latc, ln2e, ntrac, late
  logical :: semice, qflux, nsib, qcloud, coupled_aero

  ! Get names of input and output files
  call getarg(1, infile)
  call getarg(2, outfile)

  ! Open input and output files
  open(unit=inunit, file=infile, status='old', form='unformatted', &
       action='read')
  open(unit=outunit, file=outfile, status='new', form='unformatted', &
       action='write')

  ! Read headers from input file
  read (inunit) exptyp
  read (inunit) ndays1, mstep
  read (inunit) header2

  ! Get new date
  write (*, *)
  write (*, '(A,I6)') "Original value of ndays : ", ndays1
  write (*, '(A)', advance='no') "Please enter new value  : "
  read (*, *) ndays2
  write (*, *)

  ! Write headers to output file
  write (outunit) exptyp
  write (outunit) ndays2, mstep
  write (outunit) header2

  ! Determine resolution from headers
  hres = header2(6:8)
  vres = header2(10:12)
  write (*, *) "Horizontal resolution = ", hres
  write (*, *) "Vertical resolution =   ", vres

  ! Determine whether restart data for various schemes is present
  if (exptyp(5:10) == "SEMICE") then
     semice = .TRUE.
     write (*, *) "Restart file contains SEMICE data."
  else
     write (*, *) "Restart file does not contain SEMICE data."
     semice = .FALSE.
  end if
  if (header2(26:30) == "OCHFA") then
     qflux = .TRUE.
     write (*, *) "Restart file contains QFLUX data."
  else
     qflux = .FALSE.
     write (*, *) "Restart file does not contain QFLUX data."
  end if
  if (header2(38:42) == "SOIL9") then
     nsib = .TRUE.
     write (*, *) "Restart file contains NSIB data."
  else
     nsib = .FALSE.
     write (*, *) "Restart file does not contain NSIB data."
  end if
  if (exptyp(47:50) == "QCL4") then
     qcloud = .TRUE.
     write (*, *) "Restart file contains QCLOUD data."
  else
     write (*, *) "Restart file does not contain QCLOUD data."
     qcloud = .FALSE.
  end if
  if (exptyp(12:14) == "AER") then
     coupled_aero = .TRUE.
     write (*, *) "Restart file contains aerosol data."
  else
     write (*, *) "Restart file does not contain aerosol data."
     coupled_aero = .FALSE.
  end if

  ! Allocate data arrays accordingly
  select case (hres)
  case ("R21")
     lon = 64
     lat = 28
     nw = 484
  case ("R42")
     lon = 128
     lat = 54
     nw = 1849
  case ("T63")
     lon = 192
     lat = 48
     nw = 2080
  case default
     stop 'Invalid horizontal resolution.'
  end select
  select case (vres)
  case (" 9L")
     nl = 9
  case ("18L")
     nl = 18
  case ("24L")
     nl = 24
  case default
     stop 'Invalid vertical resolution.'
  end select
  ln2 = 2 * lon
  lat2 = 2 * lat
  lonc = ln2
  latc = lat
  ln2e = ln2
  ntrac = 3
  late = lat
  allocate (spec(nw), qout(ln2, lat), imsl(ln2, lat), griddata(ln2, ngrid), &
            ron(ln2, nl), son(ln2, nl), sdot(lon, nl, lat2), sgsav(ln2, lat))
  if (semice) allocate (stice(ln2, 9))
  if (qflux) allocate (ochfa(ln2, lat))
  if (nsib) allocate (wb(ln2, ms, lat), tggsn(ln2, 3, lat), ssdnn(ln2, lat), &
                      isflag(ln2, lat))
  if (qcloud) allocate (qlb(lonc, nl, latc), oprec1(ln2, lat))
  if (coupled_aero) allocate (bigxt(ln2e, nl, ntrac, late))

  ! Read input file, writing data immediately to output file
  n = 6 * nl + 3
  do i = 1, n
     read (inunit) spec
     write (outunit) spec
  end do

  n = 3 * nl
  do i = 1, n
     read (inunit) qout
     write (outunit) qout
  end do

  read (inunit) imsl
  write (outunit) imsl

  do i = 1, lat
     read (inunit) griddata
     write (outunit) griddata
     read (inunit) ron, son
     write (outunit) ron, son
  end do

  if (semice) then
     do i = 1, lat
        read (inunit) stice
        write (outunit) stice
     end do
  end if

  n = 2 * nl
  do i = 1, n
     read (inunit) spec
     write (outunit) spec
  end do

  read (inunit) sdot
  write (outunit) sdot

  if (qflux) then
     read (inunit) ochfa
     write (outunit) ochfa
  end if

  if (nsib) then
     do i = 1, 3
        read (inunit) wb
        write (outunit) wb
     end do
     read (inunit) tggsn
     write (outunit) tggsn
     read (inunit) ssdnn
     write (outunit) ssdnn
     do i = 1, 2
        read (inunit) tggsn
        write (outunit) tggsn
     end do
     do i = 1, 2
        read (inunit) ssdnn
        write (outunit) ssdnn
     end do
     read (inunit) isflag
     write (outunit) isflag
     do i = 1, 8
        read (inunit) ssdnn
        write (outunit) ssdnn
     end do
  end if

  do i = 1, 6
     read (inunit) sgsav
     write (outunit) sgsav
  end do

  if (qcloud) then
     do i = 1, 5
        read (inunit) qlb
        write (outunit) qlb
     end do
     read (inunit) oprec1
     write (outunit) oprec1
  end if

  if (coupled_aero) then
     do i = 1, 2
        read (inunit) bigxt
        write (outunit) bigxt
     end do
  end if


  ! Close input and output files
  close (unit=inunit)
  close (unit=outunit)

end program redate_restart_mk3