PROGRAM rejouer

  USE mod_const_mpi, ONLY: comm_lmdz
  USE inigeomphy_mod, ONLY: inigeomphy
  USE comvert_mod, ONLY: presnivs
  USE comvert_mod, ONLY: preff, pa
  USE ioipsl, ONLY: getin

  IMPLICIT NONE
  INCLUDE "dimensions.h"

  REAL :: airefi
  REAL :: zcufi = 1.
  REAL :: zcvfi = 1.
  REAL :: rlat_rad(1), rlon_rad(1)

  INTEGER ntime
  INTEGER jour0, mois0, an0, day_step, anneeref, dayref
  INTEGER klev, klon
  CHARACTER (len = 10) :: calend
  CHARACTER(len = 20) :: calendrier


  !---------------------------------------------------------------------
  ! L'appel a inigeomphy n'est utile que pour avoir getin_p dans
  ! les initialisations
  !---------------------------------------------------------------------
  zcufi = 1.
  zcvfi = 1.
  rlat_rad(1) = 0.
  rlon_rad(1) = 0.

  preff = 101325.
  !preff=100000.
  pa = 50000.
  CALL disvert()
  CALL inigeomphy(1, 1, llm, &
          1, comm_lmdz, &
          (/rlat_rad(1), 0./), (/0./), &
          (/0., 0./), (/rlon_rad(1), 0./), &
          (/ (/airefi, 0./), (/0., 0./) /), &
          (/zcufi, 0., 0., 0./), &
          (/zcvfi, 0./))

  CALL suphel
  !ntime=4320
  ntime = 10000000
  dayref = 1
  anneeref = 2000
  CALL getin('dayref', dayref)
  CALL getin('anneeref', anneeref)
  CALL getin('calend', calend)
  CALL getin('day_step', day_step)
  calendrier = calend
  IF (calendrier == "earth_360d") calendrier = "360_day"

  jour0 = dayref
  mois0 = (jour0 - 1) / 30 + 1
  jour0 = jour0 - 30 * ((jour0 - 1) / 30)
  an0 = anneeref

  !PRINT*,"REPLAY1D jour0,mois0,an0",jour0,mois0,an0

  klon = 1
  klev = llm
  CALL iotd_ini('phys.nc', 1, 1, klev, 0., 0., presnivs, jour0, mois0, an0, 0., 86400. / day_step, calendrier)
  ! Consistent with ... CALL iophys_ini(600.)

  !---------------------------------------------------------------------
  ! Initialisation de la parametrisation
  !---------------------------------------------------------------------
  CALL call_ini_replay

  !---------------------------------------------------------------------
  ! Boucle en temps sur l'appel à la parametrisation
  !---------------------------------------------------------------------
  CALL call_param_replay(klon, klev)

end

!---------------------------------------------------------------------
!/////////////////////////////////////////////////////////////////////
!/////////////////////////////////////////////////////////////////////
!     routine additionnelles utiles, prises dans 1DUTILS.h
!/////////////////////////////////////////////////////////////////////
!/////////////////////////////////////////////////////////////////////
!---------------------------------------------------------------------

!=======================================================================
SUBROUTINE abort_gcm(modname, message, ierr)
  USE IOIPSL
  ! Stops the simulation cleanly, closing files and printing various
  ! comments
  !=======================================================================

  !  Input: modname = name of calling program
  !         message = stuff to print
  !         ierr    = severity of situation ( = 0 normal )

  CHARACTER(LEN = *) modname
  INTEGER ierr
  CHARACTER(LEN = *) message

  WRITE(*, *) 'in abort_gcm'
  CALL histclo
  WRITE(*, *) 'out of histclo'
  WRITE(*, *) 'Stopping in ', modname
  WRITE(*, *) 'Reason = ', message
  CALL getin_dump

  IF (ierr == 0) THEN
    WRITE(*, *) 'Everything is cool'
  else
    WRITE(*, *) 'Houston, we have a problem ', ierr
  endif
  STOP
END

!=======================================================================
SUBROUTINE gr_dyn_fi(nfield, im, jm, ngrid, pdyn, pfi)
  USE lmdz_ssum_scopy, ONLY: scopy

  IMPLICIT NONE
  !   passage d'un champ de la grille scalaire a la grille physique
  !=======================================================================

  !-----------------------------------------------------------------------
  !   declarations:
  !   -------------

  INTEGER im, jm, ngrid, nfield
  REAL pdyn(im, jm, nfield)
  REAL pfi(ngrid, nfield)

  INTEGER j, ifield, ig

  !-----------------------------------------------------------------------
  !   calcul:
  !   -------

  IF(ngrid/=2 + (jm - 2) * (im - 1).AND.ngrid/=1)                          &
          STOP 'probleme de dim'
  !   traitement des poles
  CALL SCOPY(nfield, pdyn, im * jm, pfi, ngrid)
  CALL SCOPY(nfield, pdyn(1, jm, 1), im * jm, pfi(ngrid, 1), ngrid)

  !   traitement des point normaux
  DO ifield = 1, nfield
    DO j = 2, jm - 1
      ig = 2 + (j - 2) * (im - 1)
      CALL SCOPY(im - 1, pdyn(1, j, ifield), 1, pfi(ig, ifield), 1)
    ENDDO
  ENDDO

  RETURN
END