;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME: read_ncdf
;
; PURPOSE:fonction de lecture pour fichier net_cdf.
; Ce programme, est moins universel que ncdf_lec (il fait appelle au
; variables declarees dans common.pro) mais il est du cop bcp plus
; facile d''utilisation. Il prend en compte la declaration des
; differents zoom qui ont ete definis (ixminmesh...premierx...) la
; declaration de la variable key_shift... bref le resultat de
; read_ncdf peut dorectement etre utilise dans plt...
; C''est aussi ce programme qui est utilise par defaut dans mes
; widgets pour la partie lecture.
;
; CATEGORY:lecture de fichiers NetCdf
;
; CALLING SEQUENCE:res = read_ncdf(name,debut[,fin])
; 
; INPUTS: name: un string definissant le champ a lire.
;         debut et fin: sont relatifs a l''axe des temps. Ce peut etre
;         - 2 dates du type yyyymmdd et ds ce cas on selectionne les
;         dates qui sont comprisent entre ces 2 dates.
;         - 2 indices qui definissent entre quel et quel pas de temps
;           on doit extraire la dimension temporelle.
;         exp: ne sert a rien!
;
; KEYWORD PARAMETERS: utilisables hors du contexte des widgets
;
;        BOITE: contient la boite sur laquelle on doit faire la
;        lecture
;        FILENAME: string contennant le nom du fichier
;        IODIRECTORY: string contennant le nom du repertoire ds lequel
;        se trouve le fichier
;        TIMESTEP:activer pour specifier que debut et fin font
;        reference a des indices de l''axe du temps et non pas a des
;        dates. 
;        TOUT: activer si on veut lire le ficher sur l''ensemble du
;        domaine sans tenir compte du sous domaine definit par boite
;        ou lon1,lon2,lat1,lat2,prof1,prof2.
;        NOSTRUCT: activer si on ne veut pas que read_ncdf retourne
;        une structure mais uniquement le tableau se rapportant au
;        champ. 
;
; OUTPUTS:une stucture lisible par litchamp.pro ou un simple tableau
; si /NOSTRUCT est active 
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:le champ doit avoir une dimension temporelle
;
; EXAMPLE:
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      15/10/1999
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION read_ncdf,name,debut,fin, pour_etre_compatible, BOITE=boite, FILENAME = filename $
             , IODIRECTORY = iodirectory, PARENTIN = parentin, TIMESTEP = timestep $
             , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, _EXTRA = ex
@common
;------------------------------------------------------------
; we find the filename.
;------------------------------------------------------------
; is parent a valid widget ?
   if keyword_set(parentin) then BEGIN
      parent = long(parentin)
      parent = parent*widget_info(parent,/managed)
   ENDIF
   if n_elements(iodirectory) EQ 0 AND n_elements(iodir) NE 0 then iodirectory = iodir
   filename = isafile(filename = filename, IODIRECTORY = iodirectory)
;------------------------------------------------------------
; ouverture du fichier nom
;------------------------------------------------------------
   if size(filename, /type) NE 7 then $
    return, report('read_ncdf cancelled')
   cdfid=ncdf_open(filename)
   contient=ncdf_inquire(cdfid)
;------------------------------------------------------------
; we check if the variable name exists in the file.
;------------------------------------------------------------
   if ncdf_varid(cdfid,name) EQ -1 then $
    return, report('variable '+name+' !C not found in the file '+filename)
;
   varcontient=ncdf_varinq(cdfid,name)
;------------------------------------------------------------
; check the time axis and the debut and fin dates
;------------------------------------------------------------
   if n_elements(debut) EQ 0 then begin
      debut = 0
      timestep = 1
   endif
   if keyword_set(timestep) then begin
      premiertps = debut
      if n_elements(fin) NE 0 then derniertps = fin ELSE derniertps = premiertps
      jpt = derniertps-premiertps+1
      time = lindgen(jpt)
   ENDIF ELSE BEGIN
      date1 = juldate(debut, /vraidate)
      if n_elements(fin) NE 0 then date2 = juldate(fin, /vraidate) ELSE date2 = date1
      if keyword_set(parent) then BEGIN
         widget_control, parent, get_uvalue = top_uvalue
         filelist = extractatt(top_uvalue,'filelist')
         currentfile = (where(filelist EQ filename))[0]
         time = (*(extractatt(top_uvalue, 'fileparameters'))[currentfile]).time_counter
         premiertps = where(time EQ date1) & premiertps = premiertps[0]
         derniertps = where(time EQ date2) & derniertps = derniertps[0]
      ENDIF ELSE BEGIN
; we find the infinite dimension
         timedim = contient.recdim
         if timedim EQ -1 then $
          return, report('the file '+filename+' as no infinite dimension. !C Use the TIMESTEP keyword')
; we find the FIRST time axis      
         timeid = 0
         repeat BEGIN           ; tant que l''on a pas trouve une variable qui n''a qu''
                                ; une dimension: la dimension infinie
            timecontient=ncdf_varinq(cdfid,timeid) ; que contient la variable
            timeid = timeid+1
         endrep until (n_elements(timecontient.dim) EQ 1 AND timecontient.dim[0] EQ contient.recdim) $
          OR timeid EQ contient.nvars+1
;
         if timeid EQ contient.nvars+1 then $
          return, report('the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword')
; we must found the time origin of the julian calendar used in the
; time axis. 
; does the attribut units exist for the variable time axis?
         timeid = timeid-1
         if timecontient.natts EQ 0 then $
          return, report('the variable '+timecontient.name+' has no attribut.!C Use the TIMESTEP keyword or add the attribut ''units'' to the variable')
         attnames = strarr(timecontient.natts)
         for attiq=0,timecontient.natts-1 do attnames[attiq]=ncdf_attname(cdfid,timeid,attiq)
         if (where(attnames EQ 'units'))[0] EQ -1 then $
          return, report('Attribut ''units'' not found for the variable '+timecontient.name+' !C Use the TIMESTEP keyword')
         ncdf_attget,cdfid,timeid,'units',value
;
; time_counter:units = "seconds since 0001-01-01 00:00:00" ;
; time_counter:units = "hours since 0001-01-01 00:00:00" ;
; time_counter:units = "days since 1979-01-01 00:00:00" ;
; time_counter:units = "months since 1979-01-01 00:00:00" ;
; time_counter:units = "years since 1979-01-01 00:00:00" ;
;
; we decript the "units" attribut to find the time origin
         value = strtrim(strcompress(string(value)), 2)
         mots = str_sep(value, ' ')
         unite = mots[0]
         depart = str_sep(mots[2], '-')
         ncdf_varget, cdfid, timeid, time
         time = long(time)
         case strlowcase(unite) of
            'seconds':time = julday(depart[1], depart[2], depart[0])+time/(long(24)*3600)
            'hours':time = julday(depart[1], depart[2], depart[0])+time/long(24)
            'days':time = julday(depart[1], depart[2], depart[0])+time
            'months':BEGIN 
               for t = 0, n_elements(time)-1  do begin
                  time[t] = julday(depart[1]+time[t], depart[2], depart[0])
               endfor
            END
            'years':BEGIN
               for t = 0, n_elements(time)-1  do begin
                  time[t] = julday(depart[1], depart[2], depart[0]+time[t])
               endfor
            END
            ELSE:return, report('The "units" attribu of the time axis must be something like: !C "seconds since 0001-01-01 ..." !C "days since 1979-01-01 ..." !C "months since 1979-01-01 ..." !C "years since 1979-01-01 ..."')
         ENDCASE
         premiertps = where(time GE date1) & premiertps = premiertps[0]
         if premiertps EQ -1 THEN $
          return, report('date 1: '+strtrim(date1, 1)+' is not found in the time axis.')
         derniertps = where(time LE date2)
         if derniertps[0] EQ -1 THEN $
          return, report('the time axis as no date before date 2: '+strtrim(date2, 1))
         derniertps = derniertps[n_elements(derniertps)-1]
      ENDELSE
      time = time[premiertps:derniertps]
      jpt = derniertps-premiertps+1
   ENDELSE
;------------------------------------------------------------
; nom de la grille a laquelle se rapporte le champ
;------------------------------------------------------------
   vargrid = strupcase(strmid(filename, strlen(filename)-9))
   if vargrid EQ 'GRID.T.NC' OR vargrid EQ 'GRID.U.NC' $
    OR vargrid EQ 'GRID.V.NC' OR vargrid EQ 'GRID.F.NC' $
    OR vargrid eq 'GRID.W.NC' $
    OR vargrid EQ 'GRID_T.NC' OR vargrid EQ 'GRID_U.NC' $
    OR vargrid EQ 'GRID_V.NC' OR vargrid EQ 'GRID_F.NC' $
    OR vargrid eq 'GRID_W.NC' then $
    vargrid = strupcase(strmid(filename, strlen(filename)-4, 1)) ELSE vargrid='T'
;---------------------------------------------------------------
; redefinition du domaine 
;---------------------------------------------------------------
   if keyword_set(tout) then begin
      nx=jpi
      ny=jpj
      nz=jpk
      premierx = 0
      premiery = 0
      premierz = 0
      dernierx = jpi-1
      derniery = jpj-1
      dernierz = jpk-1
      case strupcase(vargrid) of
         'T':mask = tmask
         'U':mask = umask()
         'V':mask = vmask()
         'W':mask = tmask
         'F':mask = fmask()
      endcase
   ENDIF ELSE BEGIN
      if keyword_set(boite) then BEGIN 
         oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
         Case 1 Of
            N_Elements(Boite) Eq 1:bte=[lon1, lon2, lat1, lat2, 0.,boite[0]]
            N_Elements(Boite) Eq 2:bte=[lon1, lon2, lat1, lat2, boite[0],boite[1]]
            N_Elements(Boite) Eq 4:bte=[Boite, prof1, prof2]
            N_Elements(Boite) Eq 5:bte=[Boite[0:3], 0, Boite[4]]
            N_Elements(Boite) Eq 6:bte=Boite
            Else: return, report('Mauvaise Definition de Boite')
         ENDCASE
         domdef, bte,GRILLE=['T', vargrid], _extra = ex
      ENDIF
      grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz
      undefine, glam & undefine, gphi & ; on libere un peu de memoire!
   ENDELSE
;---------------------------------------------------------------------
; on initialise les ixmindta, iymindta au besoin
;---------------------------------------------------------------------
   if n_elements(jpidta) EQ 0 THEN jpidta = jpiglo
   if n_elements(jpjdta) EQ 0 THEN jpjdta = jpjglo
   if n_elements(jpkdta) EQ 0 THEN jpkdta = jpkglo
   if n_elements(ixmindta) EQ 0 THEN ixmindta = 0
   if n_elements(ixmaxdta) EQ 0 then ixmaxdta = jpidta-1
   if ixmindta EQ -1 THEN ixmindta = 0
   IF ixmaxdta EQ -1 then ixmaxdta = jpidta-1
   if n_elements(iymindta) EQ 0 THEN iymindta = 0
   IF n_elements(iymaxdta) EQ 0 then iymaxdta = jpjdta-1
   if iymindta EQ -1 THEN iymindta = 0
   IF iymaxdta EQ -1 then iymaxdta = jpjdta-1
   if n_elements(izmindta) EQ 0 THEN izmindta = 0
   IF n_elements(izmaxdta) EQ 0 then izmaxdta = jpkdta-1
   if izmindta EQ -1 THEN izmindta = 0
   IF izmaxdta EQ -1 then izmaxdta = jpkdta-1
;----------------------------------------------------------------
; on va lire le fichier
;---------------------------------------------------------------
   if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1]
   key_stride = 1l > long(key_stride)
   jpitotal = long(ixmaxmesh-ixminmesh+1)
   key_shift = long(testvar(var = key_shift))
;
   key = long(key_shift MOD jpitotal)
   if key LT 0 then key = key+jpitotal
   ixmin = ixminmesh-ixmindta
;
@read_ncdf_varget
;---------------------------------------------------------------------
; on definit les variables globales rattachees a la variable
;---------------------------------------------------------------------
; varname
   varname = name
; varunit
   if varcontient.natts NE 0 then begin
      attnames = strarr(varcontient.natts)
      for attiq=0,varcontient.natts-1 do attnames[attiq]=strlowcase(ncdf_attname(cdfid,name,attiq))
;
      if (where(attnames EQ 'units'))[0] EQ -1 then value = '' $
      ELSE ncdf_attget,cdfid,name,'units',value
      varunit = strtrim(string(value), 2)
;
      if (where(attnames EQ 'add_offset'))[0] EQ -1 then add_offset = 0. $
      ELSE ncdf_attget,cdfid,name,'add_offset',add_offset
;
      if (where(attnames EQ 'scale_factor'))[0] EQ -1 then scale_factor = 1. $
      ELSE ncdf_attget,cdfid,name,'scale_factor',scale_factor
;
      if (where(attnames EQ 'missing_value'))[0] EQ -1 then missing_value = 'no' $
      ELSE ncdf_attget,cdfid,name,'missing_value',missing_value
;
   ENDIF ELSE BEGIN 
      varunit = ''
      add_offset = 0.
      scale_factor = 1.
      missing_value = 'no'
   ENDELSE
; vardate
; on construit une belle date lisible en fonction du langage specifie !!!
   nothing = juldate(debut,/vrai) ; bidouille pour recuperer year, month, day
   if keyword_set(language) then BEGIN 
      if language eq 'gb' THEN $ 
       vardate = string(format='(C(CMoA))',31*(month-1))+' '+strtrim(day, 1)+' '+strtrim(year, 1) $
      ELSE vardate = strtrim(day, 1)+' '+string(format='(C(CMoA))',31*(month-1))+' '+strtrim(year, 1)
   ENDIF ELSE vardate = strtrim(day, 1)+' '+string(format='(C(CMoA))',31*(month-1))+' '+strtrim(year, 1)
; varexp ???
; on applique la valeur valmask sur les points terre
   if NOT keyword_set(cont_nofill) then begin
      valmask = 1e20
      case 1 of
         varcontient.ndims eq 2:BEGIN ;xy array
            mask = mask[*, *, 0]
            earth = where(mask EQ 0)
            if earth[0] NE -1 then res[earth] = 1e20
         END
         varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array
            earth = where(mask EQ 0)
            if earth[0] NE -1 then res[earth] = 1e20
         END
         varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array
            mask = mask[*, *, 0]
            earth = where(mask EQ 0)
            if earth[0] NE -1 then BEGIN
               earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt)
               res[earth] = 1e20
            END
         END
         varcontient.ndims eq 4:BEGIN ;xyzt array
            earth = where(mask EQ 0)
            if earth[0] NE -1 then BEGIN
               earth = earth#replicate(1, jpt)+replicate(nx*ny*nz, n_elements(earth))#lindgen(jpt)
               res[earth] = 1e20
            END
         END
      endcase
   endif
; on applique les add_offset, scale_factor et missing_value
   if size(missing_value, /type) NE 7 then BEGIN
      if missing_value NE valmask then begin
         if abs(missing_value) LT 1e6 then missing = where(res EQ missing_value) $
         ELSE missing = where(abs(res) gt abs(missing_value)/10.)
         if missing[0] NE -1 then res[temporary(missing)] = !values.f_nan
      endif
   ENDIF
   if scale_factor NE 1 then res = res*scale_factor
   if add_offset NE 0 then res = res+add_offset
;
   if keyword_set(key_reverse) then res = reverse(temporary(res),  2)
;---------------------------------------------------------------------
   ncdf_close,cdfid
;---------------------------------------------------------------------
   if keyword_set(oldboite) then domdef, oldboite,GRILLE=['T', vargrid]
   if keyword_set(nostruct) then return, res $
   ELSE return, {tab:res, grille:vargrid, unite:varunit, experience:varexp, nom:varname}

end