Lire NetCDF.F90

!*****ISMER-UQAR********************************************************
! Code du modele ROM en developpement a ISMER-UQAR.
! POLR
!
! Revision:            $Id$
!***********************************************************************
     subroutine lire_netcdf(ncid,varid,nt,nj,ni,nk,ki,kf,datanc) 
!
!               Modification: Simon Senneville - ISMER - 5 octobre 2010 
!               Auteur: Simon Senneville - ISMER
!
! Objet:  Lecture des fichiers NetCDF
!
! Parametres:
!       intent(in):  ncid     ! Identificateur du fichier NetCDF
!       intent(in):  varid    ! Identificateur de la variable voulue
!       intent(in):  nt       ! dimension relative au temps
!       intent(in):  nj       ! dimension relative a l'axe verticale
!       intent(in):  ni       ! dimension relative a l'axe i
!       intent(in):  nk       ! dimension relative a l'axe k
!       intent(in):  ki       ! position initiale de sur l'axe KSTEP
!       intent(in):  kf       ! position finale de sur l'axe KSTEP
!       intent(out): datanc   ! champs lu
!
! Retour:       n/a              
!       [ex:    0 :  succes]
!       [      -1 :  erreur]
!
! [note: ]
!                            ::                            !
!***********************************************************************
      implicit none
! Data ditionary: declare variables type, definition and units
!     type                   :: var_name                   ! definition, units
!    [type,dimension([x,y])  :: var_name                   ! definition, units]
!    [type,intent(in|out|inout]):: name                    ! definition, units] ...
     integer,intent(in)     :: nt, ni, nj, nk, ki, kf, varid, ncid
     real,dimension(nt,nj,ni,nk),intent(out) :: datanc
     real,dimension(nk,ni)  :: datanc_no_comp
     integer                :: i, k
!-----------------------------------------------------------------------
     integer                :: idims, it, iw
     integer                :: ret, xtype, ndims, natts
     integer                :: id, wi2d, wi3d, dimw, iwid, indx, indy, indz, indice
     integer                :: dimx, dimy, dimz, dimt, stat, dimxy
     integer                :: xid, yid, zid, tid, tvid, ksvid, dtid, dimtid, dimwid
     integer                :: tread_ti, tread_tf
     integer                :: deltat_nc, dt_in,kstart, kstart_p1, dkstep
     integer,dimension(2)   :: indstart_cp, indcount_cp
     integer,dimension(6)   :: ncref                      ! Date de reference (ss,mm,hh,j,m,a)
     real(kind=8)           :: nct, ncstart, ncend, ncvalid
     real,allocatable,dimension(:,:)  :: data_wet
     integer,allocatable,dimension(:) :: dimids
     integer,allocatable,dimension(:) :: wet
     character(len=80)      :: nom 
     character(len=80),allocatable,dimension(:) :: dimnames
     logical                :: v2d, v3d, comp, vtd
                                                          ! Including the NetCDF library for fortran
#include "netcdf.inc" 
     
     vtd  = .FALSE.                                       ! Initialisation de la semaphore vtd
     ret = nf_inq_varndims(ncid,varid,ndims)              ! Trouver le nombre de dimension de la variable TAG
     if (ret.ne.NF_NOERR) then
        write(*,*) nf_strerror(ret)
        stop
     end if
     allocate (dimids(ndims),dimnames(ndims),stat=stat)   ! Allouer la dimension de dimids et dimnames
     if (stat /= 0) then
        write(*,*)'Erreur: allocate'
        stop
     end if
     ret = nf_inq_var(ncid,varid,nom,xtype,ndims,dimids(1),natts) ! Lire les informaitons relatives a la variable varid
     if (ret.ne.NF_NOERR) then
        write(*,*) nf_strerror(ret)
        stop
     end if
     do idims = 1,ndims
       ret = nf_inq_dimname(ncid,dimids(idims),nom)       ! Lire de nom des dimensions
       if (ret.ne.NF_NOERR) then
          write(*,*) nf_strerror(ret)
          stop
       end if
       dimnames(idims) = nom(:)
       if (trim(dimnames(idims))=='NVALID2D') then        ! Verifier si c'est une variable 2D compressee
         comp = .TRUE.
         v2d  = .TRUE.
         ret = nf_inq_dimid(ncid,'NVALID2D',dimwid)       ! Si oui, lire l'id de NVALID2D
         if (ret.ne.NF_NOERR) then
            write(*,*) nf_strerror(ret)
            stop
         end if
         ret = nf_inq_dimlen(ncid,dimwid,dimw)            ! lire la dimension de NVALID2D
         if (ret.ne.NF_NOERR) then
            write(*,*) nf_strerror(ret)
            stop
         end if
         ret = nf_inq_varid(ncid,'INDEXVALID2D',iwid)     ! lire l'id de INDEXVALID2D
         if (ret.ne.NF_NOERR) then
            write(*,*) nf_strerror(ret)
            stop
         end if
       end if
       if (trim(dimnames(idims))=='NVALID3D') then        ! Verifier si c'est une variable 3D compressee
         comp = .TRUE.
         v3d  = .TRUE.
         ret = nf_inq_dimid(ncid,'NVALID3D',dimwid)       ! Si oui, lire l'id de NVALID3D
         if (ret.ne.NF_NOERR) then
            write(*,*) nf_strerror(ret)
            stop
         end if
         ret = nf_inq_dimlen(ncid,dimwid,dimw)            ! lire la dimension de NVALID3D
         if (ret.ne.NF_NOERR) then
            write(*,*) nf_strerror(ret)
            stop
         end if
         ret = nf_inq_varid(ncid,'INDEXVALID3D',iwid)     ! lire l'id de INDEXVALID3D
         if (ret.ne.NF_NOERR) then
            write(*,*) nf_strerror(ret)
            stop
         end if
       end if
       if (trim(dimnames(idims))=='TROM') then            ! Verifier si c'est une variable fonction de TROM
         vtd  = .TRUE.
         ret = nf_inq_dimid(ncid,'TROM',tid)              ! Si oui, lire l'id de TROM
         if (ret.ne.NF_NOERR) then
           write(*,*) nf_strerror(ret)
           stop
         end if
         ret = nf_inq_varid(ncid,'TROM',dimt)             ! lire l'id variable TROM 
         if (ret.ne.NF_NOERR) then
           write(*,*) nf_strerror(ret)
           stop
         end if
       end if
       if (trim(dimnames(idims))=='KROM') then            ! Non-compressee
         comp  = .FALSE.
       end if
       if (trim(dimnames(idims))=='IROM') then            ! Non-compressee
         comp  = .FALSE.
       end if
       if (trim(dimnames(idims))=='ZROM') then            ! Non-compressee
         comp  = .FALSE.
       end if
     end do
     if (comp) then                                       ! Si la variable est compressee, il faut la decompresse
       allocate (wet(dimw),stat=stat)                     ! Alloue l'espace a l'axe de comrpession
       if (stat/=0) then
         write(*,*) 'Error with allocate of the variable: wet'
         stop
       end if
       allocate (data_wet(dimw,nt),stat=stat)             ! Alloue l'espace a la variable compressee
       if (stat/=0) then
         write(*,*) 'Error with allocate of the variable: data_wet'
         stop
       end if
       ret = nf_get_var_int(ncid,iwid,wet)                ! Lire la variable de l'axe de compression
       if (ret.ne.NF_NOERR) then
          write(*,*) nf_strerror(ret)
          stop
       end if
                                                          !------------------------------
       if (vtd) then                                      ! Si la variable est en fonction de TROM
         indstart_cp(1) = 1                               ! Lire une ecriture
         indstart_cp(2) = ki                              ! Le ki ieme element
         indcount_cp(1) = dimw                            ! de dimw elements 
         indcount_cp(2) = nt                              ! pour nt pas de temps
         ret = nf_get_vara_real(ncid,varid,indstart_cp(1),indcount_cp(1), data_wet(1,1)); ! Lire cette information dans data_wet
         if (ret.ne.NF_NOERR) then
            write(*,*) nf_strerror(ret)
            stop
         end if
                                                          !------------------------------
       else                                               ! Si la variable est independante du temps
         indstart_cp(1) = 1                               ! Lire une ecriture
         indstart_cp(2) = 1                               ! Le ki ieme element
         indcount_cp(1) = dimw                            ! de dimw elements 
         indcount_cp(2) = 1                               ! pour nt pas de temps
         ret = nf_get_vara_real(ncid,varid,indstart_cp(1),indcount_cp(1), data_wet(1,1)); ! Lire cette information dans data_wet
         if (ret.ne.NF_NOERR) then
            write(*,*) nf_strerror(ret)
            stop
         end if
       end if
                                                          !------------------------------
       dimxy=ni*nk;                                       ! Dimension horizontale
       do it = 1 , nt                                     ! Decompression du champs compressee
         do iw = 1 , dimw
           indice = wet(iw);
           indz   = (indice-1) / dimxy;
           indice = indice-indz*dimxy;
           indy   = (indice-1) / nk;
           indx   = indice-indy*nk;
           indy   = indy + 1 ;
           indz   = indz + 1 ;
           datanc(it,indz,indy,indx) = data_wet(iw,it);   ! La variable decompresse est lu dans datanc
         end do
       end do
       deallocate(wet,data_wet)
     else                                                 ! Si la variable ne depend pas du temps
       ret = nf_get_var_real(ncid,varid,datanc_no_comp)   ! Lire cette information dans datanc directment
                                                          ! nf_get_var_real lit les variable dans un ordre different (comme en C)
       if (ret.ne.NF_NOERR) then
          write(*,*) nf_strerror(ret)
          stop
       endif
       forall (i=1:ni,k=1:nk)                             ! Transformation de (k,i) a (i,k)
         datanc(1,1,i,k)=datanc_no_comp(k,i)              ! A_FAIRE recoder cette partie afin qu'elle soit plus robuste
       endforall
     endif
     ret = nf_close(ncid)                                 ! Fermeture du fichier ncid ouvert par getdims
     if (ret.ne.NF_NOERR) then
        write(*,*) nf_strerror(ret)
        stop
     end if
     deallocate(dimids,dimnames,stat=stat)                ! Deallocation de l'espace memoire
     if (stat /= 0) then
        write(*,*)'Dans lire_netcdf Erreur: deallocate'
        stop
     end if
                                                          !------------------------------
     end subroutine lire_netcdf