PROGRAM tostdlev

!
! Pour passer a posteriori des champs 3d du modele sur des niveaux standards
!
! L. Fairhead 2004/12
!
! 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

  implicit none

#include "netcdf.inc"

  INTEGER :: lunout, ierr, icount, ivar, nvars, nlev, len
  integer :: press_id, data_id, newdata_id
  INTEGER :: lonid, ilon, newlonid
  INTEGER :: latid, ilat, newlatid
  INTEGER :: levid, ilev, newlevid
  INTEGER :: vartype, ndims
  INTEGER :: timid, itime, newtimid
  INTEGER :: varlonid, newvarlonid
  INTEGER :: varlatid, newvarlatid
  INTEGER :: vartimid, newvartimid
  INTEGER :: varlevid, newvarlevid
  INTEGER :: varid, newvarid
  integer :: presvar_id 
  INTEGER, dimension(4) :: dimids
  integer :: natts
  CHARACTER (len=80) :: varname, attname
  INTEGER, DIMENSION(4) :: start, count
  
  REAL, DIMENSION(:), ALLOCATABLE :: lon, lat, time
  REAL, DIMENSION(:,:,:), ALLOCATABLE :: pression, champ, newchamp


! niveaux standards:
  INTEGER,parameter :: nlevstd = 17
  real       :: rlevstd(nlevstd)
  DATA rlevstd /100000., 92500., 85000., 70000.,             &
    & 60000., 50000., 40000., 30000., 25000., 20000.,        & 
    & 15000., 10000., 7000., 5000., 3000., 2000., 1000./
!
! quelques initialisations
  lunout = 6

!
! Ouverture du fichier contenant les pressions modele
  ierr = nf_open('pression.nc', NF_NOWRITE, press_id)
  IF (ierr /= NF_NOERR) then
    WRITE(lunout,*)NF_STRERROR(ierr)
    stop
  endif

!
! Ouverture du fichier contenant le champ a interpoler
  ierr = nf_open('data.nc', NF_NOWRITE, data_id)
  IF (ierr /= NF_NOERR) then
    WRITE(lunout,*)NF_STRERROR(ierr)
    stop
  ENDIF

!
! Ouverture du fichier contenant le champ interpolé
  ierr = nf_create('newdata.nc', NF_CLOBBER, newdata_id)
  IF (ierr /= NF_NOERR) then
    WRITE(lunout,*)NF_STRERROR(ierr)
    stop
  endif

!
! Début lecture fichier origine/ecriture fichier sortie
!
! Definition des dimensions:
  ierr = nf_inq_ndims(data_id,ndims)
  DO icount = 1, ndims
    ierr = nf_inq_dim(data_id, icount, varname, len)
    select case (trim(varname))
    case ('lon')
      ilon = len
!     Definition de la longitude
      ierr = nf_def_dim(newdata_id,varname,ilon, newlonid)
    case ('lat')
      ilat = len
!     Definition de la latitude
      ierr = nf_def_dim(newdata_id,varname,ilat, newlatid)
    case('presnivs')
      ilev = len
!     Definition niveaux verticaux
      ierr = nf_def_dim(newdata_id,'presnivs',nlevstd,newlevid)
    case ('time_counter')
      itime = len
!     Definition du temps
      ierr = nf_def_dim(newdata_id,varname,itime, newtimid)
    case default
      WRITE(lunout,*)'je ne reconnais pas cette dimension: ',varname
      stop
    end select
  enddo
!
! Definition des variables
!
!
!
  ierr = nf_inq_nvars(data_id, nvars)
  DO ivar = 1, nvars
    ierr = nf_inq_var(data_id, ivar, varname, vartype, ndims, dimids, natts)
    if (ierr /= 0) call handle_err(ierr)
    selectcase (trim(varname))
    case('lon')
!     definition de la longitude
      ierr = nf_inq_varid(data_id,varname, varlonid)
      dimids(1)=newlonid
      ierr = nf_def_var(newdata_id, varname, vartype , ndims, dimids,newvarlonid)
!     recopie des attributs de la variable:
      DO icount = 1, natts  
        ierr = nf_inq_attname(data_id, varlonid, icount, attname)
        ierr = nf_copy_att(data_id, varlonid, attname, newdata_id, newvarlonid) 
      enddo
    case('lat')
!     definition de la latitude
      ierr = nf_inq_varid(data_id, varname, varlatid)
      dimids(1)=newlatid
      ierr = nf_def_var(newdata_id, varname, vartype , ndims, dimids,newvarlatid)
!     recopie des attributs de la variable:
      DO icount = 1, natts  
        ierr = nf_inq_attname(data_id, varlatid, icount, attname)
        ierr = nf_copy_att(data_id, varlatid, attname, newdata_id, newvarlatid) 
      enddo
    case('time_counter')
!     definition du temps
      ierr = nf_inq_varid(data_id, varname, vartimid)
      dimids(1)=newtimid
      ierr = nf_def_var(newdata_id, varname, vartype, ndims, dimids,newvartimid)
!     recopie des attributs de la variable:
      DO icount = 1, natts  
        ierr = nf_inq_attname(data_id, vartimid, icount, attname)
        ierr = nf_copy_att(data_id, vartimid, attname, newdata_id, newvartimid) 
      enddo
    case('presnivs')
!     definition des niveaux de pression
      ierr = nf_inq_varid(data_id,varname,varlevid)
      dimids(1)=newlevid
      ierr = nf_def_var(newdata_id,  varname, vartype, ndims, dimids,newvarlevid)
!     recopie des attributs de la variable: 
      ierr = nf_copy_att(data_id, varlevid, 'units', newdata_id, newvarlevid) 
      ierr = nf_copy_att(data_id, varlevid, 'title', newdata_id, newvarlevid) 
      ierr = nf_copy_att(data_id, varlevid, 'long_name', newdata_id, newvarlevid) 
    case default
!     normalement il ne reste que la variable à interpoler 4d
      IF (ndims /= 4) then
        WRITE(lunout,*)'La variable principale du fichier n''est pas 4D'
        stop
      endif
      ierr = nf_inq_varid(data_id,varname,varid)
      if (ierr /= 0) call handle_err(ierr)
      dimids(1) = newlonid
      dimids(2) = newlatid
      dimids(3) = newlevid
      dimids(4) = newtimid
      ierr = nf_def_var(newdata_id,  varname, vartype, ndims, dimids,newvarid)
      if (ierr /= 0) call handle_err(ierr)
      DO icount = 1, natts  
        ierr = nf_inq_attname(data_id, varid, icount, attname)
        ierr = nf_copy_att(data_id, varid, attname, newdata_id, newvarid) 
      enddo
    end select
  enddo
!

!
! fermeture du mode definition du fichier
  ierr = nf_enddef(newdata_id)
!
! ecriture des variables:
! longitude
  allocate(lon(ilon))
  ierr = nf_get_var_real(data_id, varlonid, lon)
  ierr = nf_put_var_real(newdata_id, newvarlonid,lon)
!
! latitude
  allocate(lat(ilat))
  ierr = nf_get_var_real(data_id, varlatid, lat)
  ierr = nf_put_var_real(newdata_id, newvarlatid, lat)
!
! niveaux de pression
  ierr = nf_put_var_real(newdata_id, newvarlevid, rlevstd)
!
! temps
  allocate(time(itime))
  ierr = nf_get_var_real(data_id, vartimid, time)
  ierr = nf_put_var_real(newdata_id, newvartimid, time)
 
! Traitement de la variable
  ALLOCATE(pression(ilon,ilat,ilev))
  ierr = nf_inq_varid(press_id, 'pres', presvar_id)
  ALLOCATE(champ(ilon,ilat,ilev))
  ALLOCATE(newchamp(ilon,ilat,nlevstd))
  start(1) = 1
  start(2) = 1
  start(3) = 1
  COUNT(1) = ilon
  COUNT(2) = ilat
  COUNT(3) = ilev
  count(4) = 1
  DO icount = 1, itime
    COUNT(3) = ilev
    start(4) = icount
    ierr = nf_get_vara_real(press_id, presvar_id, start, count, pression)
    if (ierr /= 0) call handle_err(ierr)
    ierr = nf_get_vara_real(data_id, varid, start, count, champ)
    if (ierr /= 0) call handle_err(ierr)
    DO nlev = 1, nlevstd
      CALL plevel(ilon*ilat, ilev,.TRUE., pression, rlevstd(nlev),champ, &
 &                newchamp(:,:,nlev))
    enddo
    count(3) = nlevstd
    ierr = nf_put_vara_real(newdata_id, newvarid, start, count, newchamp)
    if (ierr /= 0) call handle_err(ierr)
  enddo

! On ferme!
  ierr = nf_close(press_id)
  ierr = nf_close(data_id)
  ierr = nf_close(newdata_id)

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 tostdlev