!
! $Header$
!
PROGRAM ts2IPCC

!
! Filtre permettant de transformer les fichiers de serie temporelle a une 
! variable de l'IPSL en fichiers acceptables par le PCMDI/IPCC. 
!
! Utilisation de la bibliothèque CMOR du PCMDI
!
! L. Fairhead 2004/08
!
! Ce programme est appelé avec un argument, le nom du fichier à traiter
! Il nécessite aussi un fichier config.def contenant diverses informations 
! (voir plus bas) et l'accès à un tableau faisant la correspondance entre 
! les noms de variables du modèle et les noms imposés par l'IPCC
! Pour l'instant on ne traite que les fichiers contenant la serie 
! temporelle d'une seule variable

  use cmor_users_functions
  use netcdf

  implicit none

#include "netcdf.inc"

  CHARACTER (len=256)        :: orig_file   ! nom du fichier à traiter
  character (len=512)        :: line_read
  CHARACTER (len=128)        :: inpath, contact, hist_gen,repert,instit 
  CHARACTER (len=128)        :: hist_var,expt_id,source,comment,refs
  CHARACTER (len=20)         :: action
  character (len=20)         :: first_part
  character (len=1004)       :: second_part
  CHARACTER (len=20), DIMENSION(100) :: ipsl_name, ipsl_units
  CHARACTER (len=20), DIMENSION(100) :: ipcc_name, ipcc_table
  CHARACTER (len=80)         :: varname, units, namedim

  INTEGER                    :: orig_file_id, nvars, ndims
  INTEGER                    :: verbos, exit_ctl, realis, indice,index_table
  INTEGER                    :: iargc, iostat, ierr
  INTEGER                    :: i,idim
  INTEGER, ALLOCATABLE, DIMENSION(:)       :: dimids,axis_ids,lendim
  INTEGER                    :: latid, lonid, vertid, timeid
  INTEGER                    :: varid, cmorvarid
  INTEGER                    :: ilat, ilon, ivert, itime
  INTEGER                    :: lunout      ! device de sortie

  logical                    :: found = .false.

  REAL, ALLOCATABLE, DIMENSION(:) :: lon, lat, vert, time
  REAL, ALLOCATABLE, DIMENSION(:) :: lon_bounds, lat_bounds
  REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: donnees
  DOUBLE PRECISION, DIMENSION(1)                :: rdate
  real                            :: missing_value

  external iargc

!
! quelques initialisations
  lunout = 6
  varname = 'xxxxxxxx'

!
! On vérifie que l'appel au programme a bien un argument:
  CALL getarg(1, orig_file)
  IF (iargc() == 0 .OR.  orig_file == '-h') then
    WRITE(lunout,*)' '
    WRITE(lunout,*)' Utilisation de ce programme: '
    WRITE(lunout,*)' ./ts2IPCC nom_de_fichier [variable]'
    WRITE(lunout,*)'        ou nom_de_fichier est le nom du fichier a traiter'
    WRITE(lunout,*)'        et variable la variable a traiter [optionel]'
    WRITE(lunout,*)' '
    WRITE(lunout,*)' ./ts2IPCC -h sort ce message'
    WRITE(lunout,*)' '
    stop
  ENDIF
  if (iargc() == 2) then
    CALL getarg(2, varname)
  endif

!
! Lecture du fichier de configuration
  OPEN (20, IOSTAT=iostat, file='config.def',form='formatted')
  IF (iostat /= 0) then
    WRITE(lunout,*)'Erreur ouverture du fichier config.def'
    stop
  endif

  do while (iostat == 0)
    READ(20,'(A)',iostat=iostat)line_read
    line_read = trim(line_read)
    IF (INDEX(line_read, '#') /= 1) THEN
      first_part = trim(line_read(1:INDEX(line_read, '=')-1))
      second_part = trim(line_read(INDEX(line_read, '=')+1:))
      selectcase(first_part)
        case('inpath')
         inpath = trim(second_part)
        case('file_action')
         action = trim(second_part)
        case('verbosity')
          READ(second_part,'(i)') verbos
        case('exit_control')
          READ(second_part,'(i)') exit_ctl
        case('repertoire')
         repert = trim(second_part)
        case('experiment_ID')
         expt_id = trim(second_part)
        case('institut')
         instit = trim(second_part)
        case('source')
         source = trim(second_part)
        case('realisation')
          READ(second_part,'(i)') realis
        case('hist_gen')
         hist_gen = trim(second_part)
        case('comment')
         comment = trim(second_part)
        case('refs')
         refs = trim(second_part)
        case('hist_var')
         hist_var = trim(second_part)
        case('contact')
          contact = trim(second_part)
      end select
    endif
  enddo
  if (iostat > 0) then
    WRITE(lunout,*)'Probleme de lecture du fichier config.def, iostat = ',iostat
    stop
  endif
  close(20)

!
! Lecture du tableau de correspondance nom IPSL <=> nom IPCC
  OPEN (20, IOSTAT=iostat, file='table.def',form='formatted')
  IF (iostat /= 0) then
    WRITE(lunout,*)'Erreur ouverture du fichier table.def'
    stop
  endif
  indice = 0
  do while (iostat == 0)
    READ(20,'(A)',iostat=iostat)line_read
    line_read = trim(line_read)
    IF (INDEX(line_read, '#') /= 1) THEN
      indice = indice + 1
      ipsl_name(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
      line_read = trim(line_read(INDEX(line_read, '|')+1:))
      ipsl_units(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
      line_read = trim(line_read(INDEX(line_read, '|')+1:))
      ipcc_name(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
      ipcc_table(indice) = trim(line_read(INDEX(line_read, '|')+1:))
    endif
  enddo
  indice = indice - 1
  close(20)
!  DO i = 1, indice
!    WRITE(lunout,*)ipsl_name(i),ipsl_units(i),ipcc_name(i),ipcc_table(i)
!  enddo

!
! Ouverture du fichier a traiter
  ierr = nf_open(orig_file, NF_NOWRITE, orig_file_id)
  IF (ierr /= NF_NOERR) then
    WRITE(lunout,*)NF_STRERROR(ierr)
    stop
  endif
!
! trouver la variable a traiter, c'est une variable a 3 ou 4 dimensions
  ierr = nf_inq_nvars(orig_file_id, nvars)
  IF (ierr /= NF_NOERR) then
    WRITE(lunout,*)NF_STRERROR(ierr)
    stop
  endif

  i = 0
  if (varname == 'xxxxxxxx') then
    DO while (.not.found)
      i = i + 1
      if (i > nvars) then
        WRITE(lunout,*)' pas de variable 3d ou 4d trouvee'
        stop
      endif
      ierr = nf_inq_varname(orig_file_id, i, varname)
      IF (ierr /= NF_NOERR) then
        WRITE(lunout,*)NF_STRERROR(ierr)
        stop
      endif
      ierr =  nf_inq_varndims(orig_file_id, i, ndims)
      IF (ierr /= NF_NOERR) then
        WRITE(lunout,*)NF_STRERROR(ierr)
        stop
      endif
      if (ndims > 2) found = .true.
    enddo
  else
    ierr = nf_inq_varid(orig_file_id, varname, varid)
      IF (ierr /= NF_NOERR) then
        WRITE(lunout,*)NF_STRERROR(ierr)
        stop
      endif
  endif
    
!
! recherche de la correspondance nom IPSL <=> nom IPCC
  found = .false.
  i = 0
  do while (.not. found)
    i = i + 1
    if (i > indice) then
      WRITE(lunout,*)'La variable ',trim(varname),' n''est pas dans le tableau de correspondance table.def'
      stop
    endif  
    IF (varname == ipsl_name(i)) THEN
      index_table = i
      found = .true.
    endif
  enddo

  WRITE(lunout,*)' found variable = ', trim(varname)
  WRITE(lunout,*)' ipcc_name = ', trim(ipcc_name(index_table))


!
! Initialisation CMOR
  ierr = cmor_setup(inpath=inpath, netcdf_file_action=action,set_verbosity=verbos,&
 &                  exit_control=exit_ctl)
  IF (ierr /= 0) then
    WRITE(lunout,*)'Probleme dans cmor_setup, ierr = ', ierr
  endif

!
! Initialisation dataset
  ierr = cmor_dataset(outpath=repert,        &
 &                    experiment_id=expt_id, &
 &                    institution=instit,    &
 &                    source=source,         &
 &                    calendar='360_day',    &
 &                    realization=realis,    &
 &                    contact=contact,       &
 &                    history=hist_gen,      &
 &                    comment=comment,       &
 &                    references=refs)
  IF (ierr /= 0) then
    WRITE(lunout,*)'Probleme dans cmor_dataset, ierr = ', ierr
  endif

!
! Definition des axes

  ierr = nf90_inq_varid(orig_file_id,TRIM(varname), varid)
  if (ierr /= 0) call handle_err(ierr)
  ierr = nf90_Inquire_Variable(orig_file_id, varid, ndims = ndims)
  if (ierr /= 0) call handle_err(ierr)
  allocate (dimids(ndims))
  allocate (axis_ids(ndims))
  allocate (lendim(ndims))
  ierr = nf90_Inquire_Variable(orig_file_id, varid, dimids = dimids)
  if (ierr /= 0) call handle_err(ierr)

  do idim = 1, ndims
    ierr = nf90_Inquire_Dimension(orig_file_id, dimids(idim), &
                          name = namedim, len = lendim(idim))
    if (ierr /= 0) call handle_err(ierr)
    units=' ' 
    selectcase(trim(namedim))
      case('lat')
!     lecture de la latitude:
        allocate(lat(lendim(idim)))
        ierr = nf_inq_varid(orig_file_id, namedim, latid)
        if (ierr /= 0) call handle_err(ierr)
        ierr = nf_get_var_real(orig_file_id, latid, lat)
        if (ierr /= 0) call handle_err(ierr)
        ierr = nf_get_att_text(orig_file_id, latid, 'units', units)
        if (ierr /= 0) call handle_err(ierr)
        allocate(lat_bounds(lendim(idim)+1))
        DO i = 2, lendim(idim)
          lat_bounds(i) = lat(i-1) - (lat(i-1) - lat(i))/2
        enddo 
        lat_bounds(1) = lat(1)
        lat_bounds(lendim(idim)+1) = lat(lendim(idim))
!       definition de la latitude
        axis_ids(idim) = cmor_axis(                           &
                    table=trim(ipcc_table(index_table)),      &
                    table_entry='latitude',                   &
                    units=units,                              &  
                    length=lendim(idim),                      &
                    coord_vals=lat,                           &
                    cell_bounds=lat_bounds)
!
!        
      case('lon')
!       lecture de la longitude:
        allocate(lon(lendim(idim)))
        ierr = nf_inq_varid(orig_file_id, namedim, lonid)
        if (ierr /= 0) call handle_err(ierr)
        ierr = nf_get_var_real(orig_file_id, lonid, lon)
        if (ierr /= 0) call handle_err(ierr)
        ierr = nf_get_att_text(orig_file_id, lonid, 'units', units)
        if (ierr /= 0) call handle_err(ierr)
        ALLOCATE(lon_bounds(lendim(idim)+1))
        DO i = 2, lendim(idim)
         lon_bounds(i) = lon(i-1) - (lon(i-1) - lon(i))/2
        enddo 
        lon_bounds(1) = lon(1) - (lon_bounds(3) -lon_bounds(2))/2.
        lon_bounds(lendim(idim)+1) = lon(lendim(idim)) + (lon_bounds(lendim(idim))-lon_bounds(lendim(idim)-1))/2. 

!       definition de la longitude
        axis_ids(idim) = cmor_axis(                           &
                    table=trim(ipcc_table(index_table)),      &
                    table_entry='longitude',                  &
                    units=units,                              &  
                    length=lendim(idim),                            &
                    coord_vals=lon,                           &
                    cell_bounds=lon_bounds)        
!
!
      case('presnivs')
!     lecture de la verticale:
        allocate(vert(lendim(idim)))
        ierr = nf_inq_varid(orig_file_id, namedim, vertid)
        if (ierr /= 0) call handle_err(ierr)
        ierr = nf_get_var_real(orig_file_id, vertid, vert)
        if (ierr /= 0) call handle_err(ierr)
        ierr = nf_get_att_text(orig_file_id, vertid, 'units', units)
        if (ierr /= 0) call handle_err(ierr)
!
!       definition de la verticale
        if (units == 'mb') units='Pa'
        axis_ids(idim) = cmor_axis(                           &
                    table=trim(ipcc_table(index_table)),      &
                    table_entry='pressure',                   &
                    units=units,                              &  
                    length=lendim(idim),                      &
                    coord_vals=vert)
!
!        
      case('time_counter')         
!     definition du temps
      if (idim /= ndims) then
        write(lunout,*)'la dimension temps doit etre la derniere dimension'
        stop
      endif
      allocate(time(lendim(idim)))
      ierr = nf_inq_varid(orig_file_id,namedim,timeid)
      if (ierr /=0) call handle_err(ierr)
      ierr = nf_get_var_real(orig_file_id, timeid, time)
      if (ierr /=0) call handle_err(ierr)
      ierr = nf_get_att_text(orig_file_id,timeid, 'units', units)
      if (ierr /=0) call handle_err(ierr)
      axis_ids(idim) = cmor_axis(                          &
           table=trim(ipcc_table(index_table)),            &
           table_entry='time',                             &
           units=units,                                    &
           length=lendim(idim),                            &
           interval='30 minutes')
      itime= lendim(idim)
    case default
      write(lunout,*)'Dimension: ', trim(namedim),' non reconnue'
      stop
   endselect
  enddo
  
!
! Definition de la variable a ecrire
  units=' ' 
  ierr = nf_inq_varid(orig_file_id,TRIM(varname), varid)
  if (ierr /= 0) call handle_err(ierr)
  ierr = nf_get_att_text(orig_file_id, varid, 'units', units) 
  if (ierr /= 0) call handle_err(ierr)
  ierr = nf_get_att_real(orig_file_id, varid, 'missing_value', missing_value)
  if (ierr /= 0) call handle_err(ierr)
  cmorvarid = cmor_variable(                         &
       table=trim(ipcc_table(index_table)),          &
       table_entry=trim(ipcc_name(index_table)),     &
       units=units,                                  &
       axis_ids= axis_ids,                           &
       missing_value=real(missing_value),            &
       original_name=varname)
!
! Lecture de la variable
  if (ndims == 3) then
    ALLOCATE (donnees(lendim(1), lendim(2), 1, lendim(3) ))
  else if (ndims ==4) then
    allocate (donnees(lendim(1), lendim(2), lendim(3), lendim(4) ))
  endif
  ierr = nf_get_var_real(orig_file_id, varid, donnees)
!
! Ecriture de la variable
 
  DO i = 1, itime
    rdate(1) = dble(time(i))
    ierr = cmor_write(                                     &
             var_id        = cmorvarid,                    &
             DATA          = REAL(donnees(:,:,:,i)),       &
             ntimes_passed = 1,                            &
             time_vals     = rdate)
  enddo
!
! Fin CMOR
  ierr = cmor_close()
  IF (ierr /= 0) then
    WRITE(lunout,*)'Probleme dans cmor_close, ierr = ', ierr
  endif

!
! fermeture fichier originel
  ierr = nf_close(orig_file_id)
  IF (ierr /= NF_NOERR) then
    WRITE(lunout,*)NF_STRERROR(ierr)
    stop
  endif

contains

subroutine handle_err(status)
  use netcdf

  implicit none
  integer, intent(in) :: status

  write(lunout,*)nf90_strerror(status)
  stop

end subroutine handle_err
  
END PROGRAM ts2IPCC