;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:ncdf_meshlec
;
; PURPOSE:lit le fichier meshmask au format netcdf cree par OPA ds sa
; version tout NETCDF!
;
; CATEGORY:lecture de grille
;
; CALLING SEQUENCE:ncdf_meshlec [,' nomfich']
; 
; INPUTS:
;
;       nomfich: un string definissant le nom du fichier a lire. Si
; ce string ne contient pas l''adresse entiere du fichier nomfich est
; complete par iodir. Par defaut nomfich=iodir+'meshmask.nc'
;
; KEYWORD PARAMETERS:
;
;    GLAMBOUNDARY: un vecteur de 2 elements specifaint le min et le
;    max qui doivent etre imposes en longitude (obligatoire si le
;    tableau depasse 360 degres)
;
;    /CHECKDAT: cherche si il existe un ficher du meme nom que le
;    fichier meshmask mais terminant par .dat au lieu de .dat. 
;             Si ce fichier n''existe pas, ncdf_meshlec lit le fichier
;         NetCdf et sauve tous les tableaux lus ds le fichier
;         ...dat. Rq: Ce fichier est bcp plus petit que le fichier
;         NetCdf d''origine car par ex on ne sauve pas les tableaux
;         umask, vmask, fmask mais simplement certains de leurs bords
;         (cf. umask.pro vmask.pro et fmask.pro) et ces masks sont
;         sauves en bytes et toues les autres tableaux en
;         float. Cependant ce fichier .dat est cree en fonction des
;         parametres rentres ds le init... (ixminmesh.....) et Je ne
;         sais rien sur la portabilite du fichier .dat.
;            Si ce fichier .dat existe on le lit avec restore.
;
; OUTPUTS:mise a jours des variables de common glam, common
; gphi,common e12,common mask,common tab (cf. common.pro)
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:prend en compte (si il sont definits) les
; ixminmesh,...et definit jpiglo,jpi,....
;
; RESTRICTIONS:
;
;  La definition de ixminmesh,ixmaxmesh,iyminmesh,iymaxmesh
;  ,izminmesh,izmaxmesh doit etre faite avant l''entree dans cette
;  routine. pour attribuer automatiquement ces valeurs au maximum
;  possible les mettre toutes a -1 et ncdf_meshlec les calculera.
;
; EXAMPLE:
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      12/1999
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
PRO ncdf_meshlec, nomfich, GLAMBOUNDARY = glamboundary, GETDIMENSIONS = getdimensions, CHECKDAT = checkdat
@common
   tempsun = systime(1)         ; pour key_performance
;----------------------------------------------------------
;       constitution de l''adresse s_fichier et
;     ouverture du fichier a l''adresse s_fichier
;-------------------------------------------------------
; def de nomfich par defaut
   IF n_params() EQ 0 then nomfich = 'meshmask.nc'
;------------------
   if keyword_set(CHECKDAT) then begin
      nomfichdat = strmid(nomfich, 0, strlen(nomfich)-2)+'dat'
      thisOS = strupcase(strmid(!version.os_family, 0, 3))
      CASE thisOS OF 
         'MAC':sep = ':'
         'WIN':sep = '\'
         ELSE:sep = '/'
      ENDCASE
      if strpos(nomfichdat, sep) EQ -1 then begin
; si iodir ne finit pas par sep on le rajoute
         if rstrpos(iodir, sep) NE strlen(iodir)-1 then iodir = iodir+sep
         nomfichdat = iodir+nomfichdat
      endif
      test = findfile(nomfichdat)  ; le nom cherche correspond bien a un fichier?
      if test[0] NE '' then begin
         noticebase=xnotice('Lecture du fichier !C '+nomfichdat+'!C ...')
         restore, nomfichdat
         widget_control, noticebase, bad_id = nothing, /destroy
         if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun 
         return
      ENDIF
   endif
;------------------
   s_fichier = isafile(file = nomfich, iodir = iodir)
;
   noticebase=xnotice('Lecture du fichier !C '+s_fichier+'!C ...')
;
   cdfid=ncdf_open(s_fichier)
   contient=ncdf_inquire(cdfid)

;------------------------------------------------------------
; differentes dimensions
;------------------------------------------------------------
   ncdf_diminq,cdfid,'x',name,jpiglo
   ncdf_diminq,cdfid,'y',name,jpjglo
   ncdf_diminq,cdfid,'z',name,jpkglo
;
   if keyword_set(getdimensions) then begin
      ncdf_close,  cdfid
      return
   endif
;-------------------------------------------------------
;-------------------------------------------------------
   if n_elements(ixminmesh) EQ 0 THEN ixminmesh = 0
   if n_elements(ixmaxmesh) EQ 0 then ixmaxmesh = jpiglo-1
   if ixminmesh EQ -1 THEN ixminmesh = 0
   IF ixmaxmesh EQ -1 then ixmaxmesh = jpiglo-1
   if n_elements(iyminmesh) EQ 0 THEN iyminmesh = 0
   IF n_elements(iymaxmesh) EQ 0 then iymaxmesh = jpjglo-1
   if iyminmesh EQ -1 THEN iyminmesh = 0
   IF iymaxmesh EQ -1 then iymaxmesh = jpjglo-1
   if n_elements(izminmesh) EQ 0 THEN izminmesh = 0
   IF n_elements(izmaxmesh) EQ 0 then izmaxmesh = jpkglo-1
   if izminmesh EQ -1 THEN izminmesh = 0
   IF izmaxmesh EQ -1 then izmaxmesh = jpkglo-1
;
   ixmindtasauve = testvar(var = ixmindta)
   iymindtasauve = testvar(var = iymindta)
   izmindtasauve = testvar(var = izmindta)
;
   ixmindta = 0l
   iymindta = 0l
   izmindta = 0l
;
   jpi    = long(ixmaxmesh-ixminmesh+1)
   jpj    = long(iymaxmesh-iyminmesh+1)
   jpk    = long(izmaxmesh-izminmesh+1)
;
   if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1]
   key_stride = 1l > long(key_stride)
;
   jpitotal = jpi
   key_shift = long(testvar(var = key_shift))
   key = long(key_shift MOD jpitotal)
   if key LT 0 then key = key+jpitotal
   ixmin = ixminmesh-ixmindta
;
   jpi = jpi/key_stride[0]
   jpj = jpj/key_stride[1]
   jpk = jpk/key_stride[2]
;
   jpt = 1
   premiertps = 0
;
   premierx = 0
   dernierx = jpi-1
   premiery = 0
   derniery = jpj-1
   premierz = 0
   dernierz = jpk-1
   nx = jpi
   ny = jpj
   nz = 1
   izminmeshsauve = izminmesh
   izminmesh = 0
;
;-------------------------------------------------------
; tableaux 2d:
;-------------------------------------------------------
   if total(key_stride) EQ 3 then $
    namevar = ['glamt','glamu','glamv','glamf','gphit','gphiu','gphiv','gphif','e1t','e1u','e1v','e1f','e2t','e2u','e2v','e2f'] $
   ELSE namevar = ['glamt','glamu','glamv','gphit','gphiu','gphiv']
;
   for i = 0, n_elements(namevar)-1 do begin
      varcontient=ncdf_varinq(cdfid, namevar[i])
      name = varcontient.name
@read_ncdf_varget
      commande = namevar[i]+'=float(res)'
      rien = execute(commande)
   ENDFOR
;
   if total(key_stride) NE 3 then BEGIN
      glamf = glamt
      glamf = [glamf, reform(2.*reform(glamt[jpi-1, *])-reform(glamt[jpi-2, *]), 1, jpj)]
      glamf = [ [glamf], [reform(2*glamf[*, jpj-1]-glamf[*, jpj-2], jpi+1, 1)] ]
      glamf = glamf+shift(glamf, -1, -1)
      glamf = .5*glamf[0:jpi-1, 0:jpj-1]
;
      gphif = gphit
      gphif = [gphif, reform(2.*reform(gphit[jpi-1, *])-reform(gphit[jpi-2, *]), 1, jpj)]
      gphif = [ [gphif], [reform(2*gphif[*, jpj-1]-gphif[*, jpj-2], jpi+1, 1)] ]
      gphif = gphif+shift(gphif, -1, -1)
      gphif = .5*gphif[0:jpi-1, 0:jpj-1]
;
      e1t = replicate(1, jpi, jpj)
      e2t = e1t
      e1u = e1t
      e2u = e1t
      e1v = e1t
      e2v = e1t
      e1f = e1t
      e2f = e1t
   ENDIF
;-------------------------------------------------------
; tableaux 3d:
;-------------------------------------------------------
   nz = jpk
   izminmesh = izminmeshsauve
;
   varcontient=ncdf_varinq(cdfid,'tmask')
   name = varcontient.name
@read_ncdf_varget
   tmask = byte(res)
;
   varcontient=ncdf_varinq(cdfid,'umask')
   name = varcontient.name
   nx = 1
   premierx = jpi-1
@read_ncdf_varget
   umaskred = byte(res)
   umaskred = reform(umaskred, /over)
;
   varcontient=ncdf_varinq(cdfid,'vmask')
   name = varcontient.name
   nx = jpi
   premierx = 0
   ny = 1
   premiery = jpj-1
@read_ncdf_varget
   vmaskred = byte(res)
   vmaskred = reform(vmaskred, /over)
;
   varcontient=ncdf_varinq(cdfid,'fmask')
   name = varcontient.name
   nx = 1
   premierx = jpi-1
   ny = jpj
   premiery = 0
@read_ncdf_varget
   fmaskredy = byte(res)
   coast = where(fmaskredy NE 0 and fmaskredy NE 1)
   IF coast[0] NE -1 THEN fmaskredy[coast] = 0b
   fmaskredy= reform(fmaskredy, /over)
   nx = jpi
   premierx = 0
   ny = 1
   premiery = jpj-1
@read_ncdf_varget
   fmaskredx = byte(res)
   coast = where(fmaskredx NE 0 and fmaskredx NE 1)
   IF coast[0] NE -1 THEN fmaskredx[coast] = 0b
   fmaskredx = reform(fmaskredx, /over)
;-------------------------------------------------------
; tableaux 1d:
;-------------------------------------------------------
   namevar = ['e3t','e3w','gdept','gdepw']
   for i = 0, n_elements(namevar)-1 do begin
      varcontient=ncdf_varinq(cdfid, namevar[i])
      if n_elements(varcontient.dim) EQ 4 THEN BEGIN
         commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $
          +',offset = [0,0,izminmesh,0], count = [1,1,jpk,1]'
         if key_stride[2] NE 1 then commande = commande+', stride=[1,1,key_stride[2],1]'
      ENDIF ELSE BEGIN
         commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $
          +',offset = [izminmesh], count = [jpk]'
         if key_stride[2] NE 1 then commande = commande+', stride=key_stride[2]'
      ENDELSE
      rien = execute(commande)
      commande = namevar[i]+'=float('+namevar[i]+')'
      rien = execute(commande)
      commande = namevar[i]+' = reform('+namevar[i]+', /over)'
      rien = execute(commande)
   ENDFOR
;-------------------------------------------------------
   ncdf_close,  cdfid
;-------------------------------------------------------
; bornes de glam qui ne doivent pas depasser 360 degres... 
;-------------------------------------------------------
;    minglam = min(glamt, max = maxglam)
;    if maxglam-minglam GE 360 AND NOT keyword_set(glamboundary) then $
;     nothing = execute('glamboundary = '+xquestion('What are the longitudes boundary?', '[-180,180]',/chkwidget))
;-------------------------------------------------------
   if keyword_set(glamboundary) then BEGIN
      if glamboundary[0] NE glamboundary[1] then BEGIN
         glamt = glamt MOD 360
         smaller = where(glamt LT glamboundary[0])
         if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360
         bigger = where(glamt GE glamboundary[1])
         if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360
         glamu = glamu MOD 360
         smaller = where(glamu LT glamboundary[0])
         if smaller[0] NE -1 then glamu[smaller] = glamu[smaller]+360
         bigger = where(glamu GE glamboundary[1])
         if bigger[0] NE -1 then glamu[bigger] = glamu[bigger]-360
         glamv = glamv MOD 360
         smaller = where(glamv LT glamboundary[0])
         if smaller[0] NE -1 then glamv[smaller] = glamv[smaller]+360
         bigger = where(glamv GE glamboundary[1])
         if bigger[0] NE -1 then glamv[bigger] = glamv[bigger]-360
         glamf = glamf MOD 360
         smaller = where(glamf LT glamboundary[0])
         if smaller[0] NE -1 then glamf[smaller] = glamf[smaller]+360
         bigger = where(glamf GE glamboundary[1])
         if bigger[0] NE -1 then glamf[bigger] = glamf[bigger]-360
      endif
   endif
;-------------------------------------------------------
   ixmindta = ixmindtasauve
   iymindta = iymindtasauve
   izmindta = izmindtasauve
;-------------------------------------------------------
   widget_control, noticebase, bad_id = nothing, /destroy
;
   if keyword_set(checkdat) then BEGIN
   noticebase=xnotice('Ecriture du fichier !C '+nomfichdat+'!C ...')
;
      save, jpi, jpj, jpk, ixminmesh, ixmaxmesh, iyminmesh, iymaxmesh, izminmesh, izmaxmesh, glamt,glamu,glamv,glamf,gphit,gphiu,gphiv,gphif,e1t,e1u,e1v,e1f,e2t,e2u,e2v,e2f, tmask, umaskred, vmaskred, fmaskredx, fmaskredy, gdept, gdepw, e3t, e3w, filename = nomfichdat
;
      widget_control, noticebase, bad_id = nothing, /destroy
   endif
;
   if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun 

;-------------------------------------------------------
   return
end