source: trunk/SRC/Obsolete/lec.pro @ 378

;+
;
; @file_comments
; lit les fichiers Vairmer en sort:
; un tableau 2d ou 3d en fonction de nomchamp qui est le nom
; du champ a extraire (2d s'il commence par SO et 3d s'il commence par VO)
; cette fonction modifie aussi les variables globales:
; varname: trois lettres: nom de l'experience
; vargrid: nom de la grille
; vardate: date (yy)yymmdd
; varexp: nom Vairmer du champ a tracer
;
; @obsolete
;
; @categories
; Graphics, Reading
;
; @examples
;
;   IDL> resultat=lec('nom_Vairmer'[,date[,'nom_experience']])
;
; @param nomchamp {in}{required}
; 2 choix possibles:
;             1) nom de champ Vairmer (chaine de 8 caracteres en majuscule ou
; minuscule commencant par vo ou so). Dans cette methode on saute directement
; d'en-tete en en-tete jusqu'a trouver le bon fichier.
;             2) chaine de caracteres commencant par vo ou so suivit du
; numero de champ a aller chercher (par ex 'vo5'). Cette methode est un peu
; plus rapide car elle va directement chercher le fichier qui nous interesse.
;
; @param date {in}{optional}
; nombres de 6 ou 8 chiffres (anneemoisjour, par ex:19980507)
;
; @param nomexp {in}{optional}
; trois lettres designant le nom de l'experience
;
; @keyword ANOM {in}
; type du fichier vairmer par rapport auquel on doit calculer
;             l'anomalie ('EX','AN','SE','MO','')
;
; @keyword ECRIT {in}
; permet d'imprimer tous les noms vairmer que contient le fichier.
; ds ce cas en input on met seulement 'vo' ou 'so' la fonction retourne le
; nombre de fichiers lus.
;
; @keyword BOITE
;
; @keyword EXPANOM {in}
; si on calcule l'anom par rapport a une exper differente
;
; @keyword FILENAME
; string pour passer directement le nom du champ sans
; utiliser les inputs: nom_Vairmer',date,'nom_experience'. Rq si
; ces inputs sont qd meme donnes ils ne sont pas modifies par
; filename.
;
; @keyword GRID
; lorsque ce mot clef est active, lec retourne la liste
; des types de grilles (T, U...) auxquelles se rapportent les
; variables. ds ce cas en input on met seulement 'vo' ou 'so'.
;
; @keyword NAME
; lorsque ce mot clef est active, lec retourne la liste
; des noms des variables. ds ce cas en input on met seulement
; 'vo' ou 'so'.
;
; @keyword TOUT
; oblige lec a lire le champ sur tout le domaine qui a
; etait selectionne pour la cession en cours (jpi,jpj,jpk)
;
; @returns
; un tableau 2d ou 3d. sans le mot cle /TOUT, sa taille est
; celle du sous domaine definit par <pro>domdef</pro>(nx,ny,nz). avec /TOUT le
; champ a la taille du  domaine qui a etait selectionne pour la
; cession en cours (jpi,jpj,jpk).
; pour les sous domaines cf:
;        <a href="http://www.ipsl.jussieu.fr/~smlod/sousdomaine.html">
; Retourne -1 en cas d'erreur.
;
; @uses
; <pro>common</pro>
; <pro>isnumber</pro>
; <pro>ficdate</pro>
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                       26/5/98
;                       Jerome Vialard : adaptation au format vairmer
;                                        keyword anom et expanom
;                       1/7/98
;                       Sebastien Masson (masque des terres)
;                       14/8/98
;                       Sebastien Masson (decoupe pour les sous domaines...)
;                       2/99
;
; @version
; $Id$
;
;-
FUNCTION lec, nomchamp,date,nomexp $
            , ECRIT=ecrit, ANOM=anom, BOITE=boite, EXPANOM=expanom $
            , TOUT=tout, GRID=grid, NAME=name, filename=FILENAME
;
  compile_opt idl2, strictarrsubs, obsolete
;
@common
   tempsun = systime(1)         ; pour key_performance
   z = -1
;
   if keyword_set(filename) then BEGIN
      CASE strupcase(strmid(!version.os_family, 0, 3)) of
         'MAC':sep = ':'
         'WIN':sep = '\'
         ELSE:sep = '/'
      ENDCASE
      fname = strmid(filename, rstrpos(filename, sep)+1)
      if n_elements(nomchamp) EQ 0 then nomchamp = strmid(fname,6, 2)
      if n_elements(date) EQ 0 then date = long(strmid(fname,8))
      if n_elements(nomexp) EQ 0 then nomexp = strmid(fname,0, 3)
   endif
;
   nomchamp=strupcase(nomchamp)
   dim=string(format='(a2)',nomchamp)
;print, 'nom de l''experience: ',nomchamp
;------------------------------------------------------------
; specification de la date et de l'experience
;------------------------------------------------------------
   case n_params() OF
      0:BEGIN
         if keyword_set(filename) then begin
            rien=juldate(date)
            prefix=nomexp
         ENDIF ELSE return, report('Donnez un argument en entree ou utilisez le mot clef FILENAME')
      END
      1:date=long(day)+long(month)*100+long(year)*10000
      2:rien=juldate(date)
      3:begin
         rien=juldate(date)
         prefix=nomexp
      end
   endcase
;------------------------------------------------------------
; verification de la dim du fichier
;------------------------------------------------------------
   if dim ne 'SO' and dim ne 'VO' then return, report('le nom du champ doit commencer par VO ou SO')
;------------------------------------------------------------
;   constitution de l'adresse ou aller chercher le fichier
;--------------------------------------------------------------
   s_fichier=ficdate(date,dim)
;--------------------------------------------------------------------
; ouverture du fichier a l'adresse s_fichier
;--------------------------------------------------------------------
   openr, numlec, s_fichier, /get_lun,ERROR=err, /swap_if_little_endian
   if err ne 0 then begin
;  print,!err_string
      return, -1
   endif
;taille en octet du fichier
   infofichier=fstat(numlec)
;---------------------------------------------------------------------
; definition de la taille du fichier a aller chercher: jpidta,jpjdta,jpkdta...
;---------------------------------------------------------------------
   if n_elements(jpidta) EQ 0 THEN BEGIN
      if n_elements(ixmindta) EQ 0 OR n_elements(ixmaxdta) EQ 0 then $
       jpidta = jpiglo else jpidta = ixmaxdta-ixmindta+1
   endif
   if n_elements(jpjdta) EQ 0 THEN BEGIN
      if n_elements(iymindta) EQ 0 OR n_elements(iymaxdta) EQ 0 then $
       jpjdta = jpjglo else jpjdta = iymaxdta-iymindta+1
   endif
   if n_elements(jpkdta) EQ 0 THEN BEGIN
      if n_elements(izmindta) EQ 0 OR n_elements(izmaxdta) EQ 0 then $
       jpkdta = jpkglo else jpkdta = izmaxdta-izmindta+1
   endif
;---------------------------------------------------------------------
; lecture des champs directement vers le champ ou l'en-tete que l'on recherche
; il faut savoir que:
;  le fortran ajoute au debut et a la fin de chaque write 4 octets de controle
;  les reels du model sont codes sur 4 octets
;  un caractere fait 1 octet
;-----------------------------------------------------------------------
;4 chaines de 8 caracteres+un tableau de reels+4 trucs de controle (pour les
; 2 write):
   if dim eq 'VO' then $
    taillebloc=4*8+long(jpidta)*jpjdta*jpkdta*4+4*4 else $
    taillebloc=4*8+long(jpidta)*jpjdta*4+4*4
;---------------------------------------------------------------------
; choix du type de lecture
;---------------------------------------------------------------------
   typelec=strmid(nomchamp,2,strlen(nomchamp))
   test=isnumber(typelec,numerochamp)
   if test eq 0 then begin
;--------------------------------------------------------------------
; 1) LECTURE DIRECTE D'EN-TETE en EN-TETE
;--------------------------------------------------------------------
      numerochamp=1
;---------------------------------------------------------------------
; lecture des noms de champ
;---------------------------------------------------------------------
      resname = ''
      resgrid = ''
      while numerochamp*taillebloc le infofichier.size do begin
         offset=(numerochamp-1)*taillebloc+4
         a=assoc(numlec,bytarr(8,/nozero), offset)
         varname=string(a[0])
         if keyword_set(ecrit) OR keyword_set(name) OR keyword_set(grid) $
          then begin
            vargrid=a[1]
            vargrid=string(vargrid[7])
            vardate=strtrim(long(string(a[2])), 2)
            varexp=strtrim(a[3], 2)
            if keyword_set(ecrit) THEN $
             print, numerochamp,' ',varname,' ',vargrid,' ',vardate,' ',varexp
            resname = [resname, varname]
            resgrid = [resgrid, vargrid]
         endif
         if nomchamp eq varname then begin
            vargrid=a[1]
            vargrid=string(vargrid[7])
            vardate=strtrim(long(string(a[2])), 2)
            varexp=strtrim(a[3], 2)
            goto,sortieboucle
         endif
         numerochamp=numerochamp+1
      ENDWHILE
      free_lun,numlec
      close, numlec
      case 1 of
         keyword_set(ecrit):return, numerochamp-1
         keyword_set(name):return, resname[1:numerochamp-1]
         keyword_set(grid):$
          return, strmid(resgrid[1:numerochamp-1],0 > (strlen(resgrid[0])-2))
         ELSE:return, report('Ce nom Vairmer de champ n''existe pas ds le fichier: '+infofichier.name)
      endcase
   endif else begin
;----------------------------------------------------------------------
; 2) LECTURE DIRECTEMENT DU CHAMP QUE L'ON VEUT
;---------------------------------------------------------------------
;---------------------------------------------------------------------
; test pour savoir si numero de champ est accessible
;---------------------------------------------------------------------
      if taillebloc*numerochamp gt infofichier.size then $
       return, report('Ce numero de champ n''existe pas. Le fichier '+infofichier.name+' ne contient que ',infofichier.size/taillebloc,' champs.')
;---------------------------------------------------------------------
; lecture de l'en-tete numero numerochamp
;---------------------------------------------------------------------
      offset=(numerochamp-1)*taillebloc+4
      a=assoc(numlec,bytarr(8,/nozero), offset)
      varname=string(a[0])
      vargrid=a[1]
      vargrid=string(vargrid[7])
      vardate=string(a[2])
      varexp=string(a[3])
   endelse
sortieboucle:
;---------------------------------------------------------------------
; lecture du champ lui-meme
;---------------------------------------------------------------------
   offset=(numerochamp-1)*taillebloc+(8+4*8)+4
   if dim eq 'VO' then $
    a=assoc(numlec,fltarr(jpidta,jpjdta,jpkdta,/nozero), offset) else $
    a=assoc(numlec,fltarr(jpidta,jpjdta,/nozero), offset)
   z=a[0]
;---------------------------------------------------------------------
; on initialise les ixmindta, iymindta  au besoin
;---------------------------------------------------------------------
   if n_elements(ixmindta) EQ 0 OR n_elements(ixmaxdta) EQ 0 then BEGIN
      ixmindta = 0
      ixmaxdta = jpidta-1
   endif
   if n_elements(iymindta) EQ 0 OR n_elements(iymaxdta) EQ 0 then BEGIN
      iymindta = 0
      iymaxdta = jpjdta-1
   endif
   if n_elements(izmin) EQ 0 OR n_elements(izmax) EQ 0 then BEGIN
      izmindta = 0
      izmaxdta = jpkdta-1
   endif
;---------------------------------------------------------------------
; on reduit z selon les valeurs de ixmindta, iymindta, ...
;---------------------------------------------------------------------
   if dim EQ 'SO' then z = z[ixminmesh-ixmindta:ixmaxmesh-ixmindta $
                             ,iyminmesh-iymindta:iymaxmesh-iymindta] $
   ELSE z = z[ixminmesh-ixmindta:ixmaxmesh-ixmindta $
              , iyminmesh-iymindta:iymaxmesh-iymindta, izminmesh-izmindta:izmaxmesh-izmindta]
;---------------------------------------------------------------------
; on shift z si key_shift est defini
;---------------------------------------------------------------------
   if n_elements(key_shift) NE 0 THEN BEGIN
      if dim EQ 'SO' then z = shift(z,key_shift, 0) $
      ELSE z = shift(z,key_shift, 0, 0)
   endif
;---------------------------------------------------------------------
;  si /TOUT n''est pas active, on coupe z pour qu''il soit a la taille
;  du zoom: nx,ny nz
;---------------------------------------------------------------------
   if NOT keyword_set(tout) then BEGIN
;-------------------------------------------------------------
; changement de domaine
;-------------------------------------------------------------
      if keyword_set(boite) then BEGIN
         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
         oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
         domdef, bte,GRILLE=vargrid
      ENDIF
;-------------------------------------------------------------
      grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz
      if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then mask = reform(mask, nx, ny, nz, /over)
      if dim EQ 'SO' then z = z[premierx:dernierx, premiery:derniery]  $
      ELSE z = z[premierx:dernierx, premiery:derniery, premierz:dernierz]
   ENDIF ELSE BEGIN
      case vargrid OF           ; on recupere le mask en entier ds le cas ou /TOUT
         'U':mask = umask()     ; n''est pas active et on le choisit en fonction
         'T':mask = tmask       ;  de la valeur de vargrid
         'W':mask = tmask
         'V':mask = vmask()
         'F':mask = fmask()
      ENDCASE
   ENDELSE
;---------------------------------------------------------------------
; calcul d'une anomalie si le keyword anom est active
;---------------------------------------------------------------------
   if keyword_set(anom) then begin
      case anom of
         'EX' : adate = 0
         'AN' : adate = floor(date/10000)*10000
         'SE' : adate = floor(date - floor(date/10000)*10000)/100 * 100
         'MO' : adate = floor(date/100)*100
         'DA' : adate = date - floor(date/10000)*10000
         ''   : adate = date - floor(date/10000)*10000
         else : return, report('Anom doit etre egal a EX,AN,SE,MO,DA ')
      endcase
      if keyword_set(expanom) then nomexpa = expanom $
      else nomexpa = nomexp
      if keyword_set(bavard) THEN print, nomchamp,' - ',adate,' - ',nomexpa
      z = z - lec(nomchamp,adate,nomexpa, TOUT = tout)
   endif
;---------------------------------------------------------------------
; on masque les terres par valmask
;---------------------------------------------------------------------
   IF n_elements(valmask) EQ 0 THEN valmask = 1e20
   if dim EQ 'SO' then BEGIN
      terre = where(mask[*,*,0] EQ 0)
      if terre[0] NE -1 then z[terre] = valmask
   ENDIF ELSE BEGIN
      terre = where(mask[*,*,0] EQ 0)
      if terre[0] NE -1 then z[where(mask EQ 0)] = valmask
   ENDELSE
;---------------------------------------------------------------------
   free_lun,numlec
   close, numlec
;---------------------------------------------------------------------
   if n_elements(oldboite) NE 0 then domdef,  oldboite
   IF keyword_set(key_performance) EQ 1 THEN print, 'temps lec', systime(1)-tempsun
;
   return,reform(z)

end