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
  
  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

  INTEGER                    :: orig_file_id, nvars, ndims
  INTEGER                    :: verbos, exit_ctl, realis, indice,index_table
  INTEGER                    :: iargc, iostat, ierr
  INTEGER                    :: i
  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

! lecture de la latitude:
  ierr = nf_inq_dimid(orig_file_id,'lat',latid)
  ierr = nf_inq_dimlen(orig_file_id, latid, ilat)
  allocate(lat(ilat))
  ierr = nf_inq_varid(orig_file_id,'lat',latid)
  ierr = nf_get_var_real(orig_file_id, latid, lat)
  ierr = nf_get_att_text(orig_file_id, latid, 'units', units)
  ALLOCATE(lat_bounds(ilat+1))
  DO i = 2, ilat
    lat_bounds(i) = lat(i-1) - (lat(i-1) - lat(i))/2
  enddo 
  lat_bounds(1) = lat(1)
  lat_bounds(ilat+1) = lat(ilat)
! definition de la latitude
  latid = cmor_axis(                         &
     table=trim(ipcc_table(index_table)),    &
     table_entry='latitude',                 &
     units=units,                            &  
     length=ilat,                            &
     coord_vals=lat,                         &
     cell_bounds=lat_bounds)        

! lecture de la longitude:
  units=' ' 
  ierr = nf_inq_dimid(orig_file_id,'lon',lonid)
  ierr = nf_inq_dimlen(orig_file_id, lonid, ilon)
  allocate(lon(ilon))
  ierr = nf_inq_varid(orig_file_id,'lon',lonid)
  ierr = nf_get_var_real(orig_file_id, lonid, lon)
  ierr = nf_get_att_text(orig_file_id, lonid, 'units', units)
  ALLOCATE(lon_bounds(ilon+1))
  DO i = 2, ilon
    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(ilon+1) = lon(ilon) + (lon_bounds(ilon)-lon_bounds(ilon-1))/2. 

! definition de la longitude
  lonid = cmor_axis(                         &
     table=trim(ipcc_table(index_table)),    &
     table_entry='longitude',                &
     units=units,                            &  
     length=ilon,                            &
     coord_vals=lon,                         &
     cell_bounds=lon_bounds)        

! definition du temps
  units=' ' 
  ierr = nf_inq_dimid(orig_file_id,'time_counter',timeid)
  ierr = nf_inq_dimlen(orig_file_id,timeid,itime)
  allocate(time(itime))
  ierr = nf_inq_varid(orig_file_id,'time_counter',timeid)
  ierr = nf_get_var_real(orig_file_id, timeid, time)
  ierr = nf_get_att_text(orig_file_id,timeid, 'units', units)
  timeid = cmor_axis(                          &
       table=trim(ipcc_table(index_table)),    &
       table_entry='time',                     &
       units=units,                            &
       length=itime,                           &
       interval='30 minutes')
  
!
! Definition de la variable a ecrire
  units=' ' 
  ierr = nf_inq_varid(orig_file_id,TRIM(varname), varid)
  ierr = nf_get_att_text(orig_file_id, varid, 'units', units) 
  ierr = nf_get_att_real(orig_file_id, varid, 'missing_value', missing_value)
  cmorvarid = cmor_variable(                         &
       table=trim(ipcc_table(index_table)),          &
       table_entry=trim(ipcc_name(index_table)),     &
       units=units,                                  &
       axis_ids=(/ lonid, latid, timeid /),          &
       missing_value=real(missing_value),            &
       original_name=varname)
!
! Lecture de la variable
  ALLOCATE (donnees(ilon, ilat, itime))
  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
  
END