;+
;
; @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 2 ou 3d. sans le mot cle /TOUT, sa taille est
; celle du sous domaine definit par domdef(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:
;
; Retourne -1 en cas d'erreur.
;
; @uses
; common
; isnumber
; fivardate
;
; @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