source: trunk/LECTURE/lec.pro @ 2

Last change on this file since 2 was 2, checked in by opalod, 22 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 15.1 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:lec
6;   
7;
8; PURPOSE:
9;       lit les fichiers Vairmer en sort:
10;               un tableau 2d ou 3d en fonction de nomchamp qui est le nom
11; du champ a extaire (2d s'il commence par SO et 3d s'il commence par VO)
12;       cette fonction modifie aussi les variables globales:
13;               varname: trois lettres: nom de l'experience
14;               vargrid: nom de la grille
15;               vardate: date (yy)yymmdd
16;               varexp: nom Vairmer du champ a tarcer
17;
18; CATEGORY:
19;       Graphics, lecture de fichier Vaimer
20;
21; CALLING SEQUENCE:
22;       resultat=lec('nom_Vairmer'[,date[,'nom_experience']])
23;
24; INPUTS:
25;       nom_Vairmer: 2 choix possibles:
26;             1) nom de champ Vairmer (chaine de 8 caracteres en majuscule ou
27; minuscule commencant par vo ou so). Dans cette methode on saute directement
28; d'en-tete en en-tete jusqu'a trouver le bon fichier.
29;             2) chaine de characteres commencant par vo ou so suivit du
30; numero de champ a aller chercher (par ex 'vo5'). Cette methode est un peu
31; plus rapide car elle va directement chercher le fichier qui nous interesse.
32;
33;       arguments optionnels:
34;       
35;       date:nombres de 6 ou 8 chiffres (anneemoisjour, par ex:19980507)
36;       nom_experience:trois lettres designant le nom de l'experience
37;
38; KEYWORD PARAMETERS:
39;
40;       /ANOM:type du fichier vairmer par rapport auquel on doit calculer
41;             l'anomalie ('EX','AN','SE','MO','')
42;
43;       /ECRIT:permet d'imprimer tous les noms vairmer que contient le fichier.
44; ds ce cas en input on met seulement 'vo' ou 'so' la fonction retourne le
45; nombre de fichiers lus.
46;
47;       /EXPANOM: si on calcule l'anom par rapport a une exper
48;       differente
49;
50;       FILENAME: string pour passer directement le nom du champ sans
51;       utiliser les inputs: nom_Vairmer',date,'nom_experience'. Rq si
52;       ces inputs sont qd meme donnes ils ne sont pas modifies par
53;       filename.
54;
55;       /GRID:lorsque ce mot clef est active, lec retourne la liste
56;       des types de grilles (T, U...) auxquelles se rapportent les
57;       variables. ds ce cas en input on met seulement 'vo' ou 'so'.
58;
59;       /NAME:lorsque ce mot clef est active, lec retourne la liste
60;       des noms des variables. ds ce cas en input on met seulement
61;       'vo' ou 'so'.
62;
63;       /TOUT: oblige lec a lire le champ sur tout le domaine qui a
64;       etait selectionne pour la cession en cours (jpi,jpj,jpk)
65;
66; OUTPUTS: un tableau 2 ou 3d. sans le mot cle /TOUT, sa taille est
67; celle du sous domaine definit par domdef (nx,ny,nz). avec /TOUT le
68; champ a la taille du  domaine qui a etait selectionne pour la
69; cession en cours (jpi,jpj,jpk).
70; pour les sous domaines cf: 
71;        http://www.ipsl.jussieu.fr/~smlod/sousdomaine.html
72;
73; COMMON BLOCKS:
74;       common.pro         isnumber.pro         fivardate.pro 
75;
76; SIDE EFFECTS:Retourne -1 en cas d'erreur.
77;
78; RESTRICTIONS:
79;
80; EXAMPLE:
81;
82; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
83;                       26/5/98
84;                       Jerome Vialard : adaptation au format vairmer
85;                                        keyword anom et expanom
86;                       1/7/98
87;                       Sebastien Masson (masque des terres)
88;                       14/8/98
89;                       Sebastien Masson (decoupe pour les sous domaines...)
90;                       2/99
91;;-
92;--------------------------------------------------------------
93;--------------------------------------------------------------
94;--------------------------------------------------------------
95function lec, nomchamp,date,nomexp,ECRIT=ecrit,ANOM=anom, BOITE = boite,EXPANOM=expanom, TOUT = tout, GRID = grid, NAME = name, filename = FILENAME
96@common
97   tempsun = systime(1)         ; pour key_performance
98   z = -1
99;
100   if keyword_set(filename) then BEGIN
101      CASE strupcase(strmid(!version.os_family, 0, 3)) of
102         'MAC':sep = ':'
103         'WIN':sep = '\'
104         ELSE:sep = '/'
105      ENDCASE
106      fname = strmid(filename, rstrpos(filename, sep)+1)
107      if n_elements(nomchamp) EQ 0 then nomchamp = strmid(fname,6, 2)
108      if n_elements(date) EQ 0 then date = long(strmid(fname,8))
109      if n_elements(nomexp) EQ 0 then nomexp = strmid(fname,0, 3)
110   endif
111;
112   nomchamp=strupcase(nomchamp)
113   dim=string(format='(a2)',nomchamp)
114;print, 'nom de l''experience: ',nomchamp
115;------------------------------------------------------------
116; specification de la date et de l'experience
117;------------------------------------------------------------
118   case n_params() OF
119      0:BEGIN
120         if keyword_set(filename) then begin
121            rien=juldate(date)
122            prefix=nomexp
123         ENDIF ELSE return, report('Donnez un argument en entree ou utilisez le mot clef FILENAME')
124      END
125      1:date=long(day)+long(month)*100+long(year)*10000
126      2:rien=juldate(date)
127      3:begin
128         rien=juldate(date)
129         prefix=nomexp
130      end
131   endcase
132;------------------------------------------------------------
133; verification de la dim du fichier
134;------------------------------------------------------------
135   if dim ne 'SO' and dim ne 'VO' then return, report('le nom du champ doit commencer par VO ou SO')
136;------------------------------------------------------------
137;   constitution de l'adresse ou aller chercher le fichier
138;--------------------------------------------------------------
139   s_fichier=ficdate(date,dim)
140;--------------------------------------------------------------------
141; ouverture du fichier a l'adresse s_fichier
142;--------------------------------------------------------------------
143   openr, numlec, s_fichier, /get_lun,ERROR=err, /swap_if_little_endian
144   if err ne 0 then begin
145;  print,!err_string
146      return, -1
147   endif
148;taille en octet du fichier
149   infofichier=fstat(numlec)
150;---------------------------------------------------------------------
151; definition de la taille du fichier a aller chercher: jpidta,jpjdta,jpkdta...
152;---------------------------------------------------------------------
153   if n_elements(jpidta) EQ 0 THEN BEGIN
154      if n_elements(ixmindta) EQ 0 OR n_elements(ixmaxdta) EQ 0 then $
155       jpidta = jpiglo else jpidta = ixmaxdta-ixmindta+1
156   endif
157   if n_elements(jpjdta) EQ 0 THEN BEGIN
158      if n_elements(iymindta) EQ 0 OR n_elements(iymaxdta) EQ 0 then $
159       jpjdta = jpjglo else jpjdta = iymaxdta-iymindta+1
160   endif
161   if n_elements(jpkdta) EQ 0 THEN BEGIN
162      if n_elements(izmindta) EQ 0 OR n_elements(izmaxdta) EQ 0 then $
163       jpkdta = jpkglo else jpkdta = izmaxdta-izmindta+1
164   endif
165;---------------------------------------------------------------------
166; lecture des champs directement vers le champ ou l'en-tete que l'on recherche
167; il faut savoir que:
168;  le fortran ajoute au debut et a la fin de chaque write 4 octets de controle
169;  les reels du model sont codes sur 4 octets
170;  un charactere fait 1 octet
171;-----------------------------------------------------------------------
172;4 chaines de 8 characteres+un tableau de reels+4 trucs de controle (pour les
173; 2 write):
174   if dim eq 'VO' then $
175    taillebloc=4*8+long(jpidta)*jpjdta*jpkdta*4+4*4 else $
176    taillebloc=4*8+long(jpidta)*jpjdta*4+4*4
177;---------------------------------------------------------------------
178; choix du type de lecture
179;---------------------------------------------------------------------
180   typelec=strmid(nomchamp,2,strlen(nomchamp))
181   test=isnumber(typelec,numerochamp)
182   if test eq 0 then begin
183;--------------------------------------------------------------------
184; 1) LECTURE DIRECTE D'EN-TETE en EN-TETE
185;--------------------------------------------------------------------
186      numerochamp=1
187;---------------------------------------------------------------------
188; lecture des noms de champ
189;---------------------------------------------------------------------
190      resname = ''
191      resgrid = ''
192      while numerochamp*taillebloc le infofichier.size do begin
193         offset=(numerochamp-1)*taillebloc+4
194         a=assoc(numlec,bytarr(8,/nozero), offset)
195         varname=string(a[0])
196         if keyword_set(ecrit) OR keyword_set(name) OR keyword_set(grid) $
197          then begin
198            vargrid=a[1]
199            vargrid=string(vargrid[7])
200            vardate=strtrim(long(string(a[2])), 2)
201            varexp=strtrim(a[3], 2)
202            if keyword_set(ecrit) THEN $
203             print, numerochamp,' ',varname,' ',vargrid,' ',vardate,' ',varexp
204            resname = [resname, varname]
205            resgrid = [resgrid, vargrid]
206         endif
207         if nomchamp eq varname then begin
208            vargrid=a[1]
209            vargrid=string(vargrid[7])
210            vardate=strtrim(long(string(a[2])), 2)
211            varexp=strtrim(a[3], 2)
212            goto,sortieboucle
213         endif
214         numerochamp=numerochamp+1
215      ENDWHILE
216      free_lun,numlec
217      close, numlec
218      case 1 of
219         keyword_set(ecrit):return, numerochamp-1
220         keyword_set(name):return, resname[1:numerochamp-1]
221         keyword_set(grid):$
222          return, strmid(resgrid[1:numerochamp-1],0 > (strlen(resgrid[0])-2))
223         ELSE:return, report('Ce nom Vairmer de champ n''existe pas ds le fichier: '+infofichier.name)
224      endcase
225   endif else begin
226;----------------------------------------------------------------------
227; 2) LECTURE DIRECTEMENT DU CHAMP QUE L'ON VEUT
228;---------------------------------------------------------------------
229;---------------------------------------------------------------------
230; test pour savoir si numero de champ est accessible
231;---------------------------------------------------------------------
232      if taillebloc*numerochamp gt infofichier.size then $
233       return, report('Ce numero de champ n''exite pas. Le fichier '+infofichier.name+' ne contient que ',infofichier.size/taillebloc,' champs.')
234;---------------------------------------------------------------------
235; lecture de l'en-tete numero numerochamp
236;---------------------------------------------------------------------
237      offset=(numerochamp-1)*taillebloc+4
238      a=assoc(numlec,bytarr(8,/nozero), offset)
239      varname=string(a[0])
240      vargrid=a[1]
241      vargrid=string(vargrid[7])
242      vardate=string(a[2])
243      varexp=string(a[3])
244   endelse
245sortieboucle:
246;---------------------------------------------------------------------
247; lecture du champ lui-meme
248;---------------------------------------------------------------------
249   offset=(numerochamp-1)*taillebloc+(8+4*8)+4
250   if dim eq 'VO' then $
251    a=assoc(numlec,fltarr(jpidta,jpjdta,jpkdta,/nozero), offset) else $
252    a=assoc(numlec,fltarr(jpidta,jpjdta,/nozero), offset)
253   z=a[0]
254;---------------------------------------------------------------------
255; on initialise les ixmindta, iymindta  au besoin
256;---------------------------------------------------------------------
257   if n_elements(ixmindta) EQ 0 OR n_elements(ixmaxdta) EQ 0 then BEGIN
258      ixmindta = 0
259      ixmaxdta = jpidta-1
260   endif
261   if n_elements(iymindta) EQ 0 OR n_elements(iymaxdta) EQ 0 then BEGIN
262      iymindta = 0
263      iymaxdta = jpjdta-1
264   endif
265   if n_elements(izmin) EQ 0 OR n_elements(izmax) EQ 0 then BEGIN
266      izmindta = 0
267      izmaxdta = jpkdta-1
268   endif
269;---------------------------------------------------------------------
270; on reduit z selon les valeurs de ixmindta, iymindta, ...
271;---------------------------------------------------------------------
272   if dim EQ 'SO' then z = z[ixminmesh-ixmindta:ixmaxmesh-ixmindta $
273                             ,iyminmesh-iymindta:iymaxmesh-iymindta] $
274   ELSE z = z[ixminmesh-ixmindta:ixmaxmesh-ixmindta $
275              , iyminmesh-iymindta:iymaxmesh-iymindta, izminmesh-izmindta:izmaxmesh-izmindta]
276;---------------------------------------------------------------------
277; on shift z si key_shift est defininit
278;---------------------------------------------------------------------
279   if n_elements(key_shift) NE 0 THEN BEGIN
280      if dim EQ 'SO' then z = shift(z,key_shift, 0) $
281      ELSE z = shift(z,key_shift, 0, 0)
282   endif
283;---------------------------------------------------------------------
284;  si /TOUT n''est pas active, on coupe z pour q''il soit a la taille
285;  du zoom: nx,ny nz
286;---------------------------------------------------------------------
287   if NOT keyword_set(tout) then BEGIN
288;-------------------------------------------------------------
289; changement de domaine
290;-------------------------------------------------------------
291      if keyword_set(boite) then BEGIN
292         Case 1 Of
293            N_Elements(Boite) Eq 1:bte=[lon1, lon2, lat1, lat2, 0.,boite[0]]
294            N_Elements(Boite) Eq 2:bte=[lon1, lon2, lat1, lat2, boite[0],boite[1]]
295            N_Elements(Boite) Eq 4:bte=[Boite, prof1, prof2]
296            N_Elements(Boite) Eq 5:bte=[Boite[0:3], 0, Boite[4]]
297            N_Elements(Boite) Eq 6:bte=Boite
298            Else: return, report('Mauvaise Definition de Boite')
299         endcase
300         oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
301         domdef, bte,GRILLE=vargrid
302      ENDIF
303;-------------------------------------------------------------
304      grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz
305      if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then mask = reform(mask, nx, ny, nz, /over)
306      if dim EQ 'SO' then z = z[premierx:dernierx, premiery:derniery]  $
307      ELSE z = z[premierx:dernierx, premiery:derniery, premierz:dernierz]
308   ENDIF ELSE BEGIN
309      case vargrid OF           ; on recupere le mask en entier ds le cas ou /TOUT
310         'U':mask = umask()     ; n''est pas active et on le choisit en fontion
311         'T':mask = tmask       ;  de la valeur de vargrid
312         'W':mask = tmask       
313         'V':mask = vmask()
314         'F':mask = fmask()
315      ENDCASE
316   ENDELSE
317;---------------------------------------------------------------------
318; calcul d'une anomalie si le keyword anom est active
319;---------------------------------------------------------------------
320   if keyword_set(anom) then begin
321      case anom of
322         'EX' : adate = 0
323         'AN' : adate = floor(date/10000)*10000
324         'SE' : adate = floor(date - floor(date/10000)*10000)/100 * 100
325         'MO' : adate = floor(date/100)*100
326         'DA' : adate = date - floor(date/10000)*10000
327         ''   : adate = date - floor(date/10000)*10000
328         else : return, report('Anom doit etre egal a EX,AN,SE,MO,DA ')
329      endcase
330      if keyword_set(expanom) then nomexpa = expanom $
331      else nomexpa = nomexp
332      if keyword_set(bavard) THEN print, nomchamp,' - ',adate,' - ',nomexpa
333      z = z - lec(nomchamp,adate,nomexpa, TOUT = tout)
334   endif
335;---------------------------------------------------------------------
336; on masque les terres par valmask
337;---------------------------------------------------------------------
338   IF n_elements(valmask) EQ 0 THEN valmask = 1e20
339   if dim EQ 'SO' then BEGIN
340      terre = where(mask[*,*,0] EQ 0)
341      if terre[0] NE -1 then z[terre] = valmask
342   ENDIF ELSE BEGIN
343      terre = where(mask[*,*,0] EQ 0)
344      if terre[0] NE -1 then z(where(mask EQ 0)) = valmask
345   ENDELSE
346;---------------------------------------------------------------------
347   free_lun,numlec
348   close, numlec
349;---------------------------------------------------------------------
350   if n_elements(oldboite) NE 0 then domdef,  oldboite
351   IF keyword_set(key_performance) EQ 1 THEN print, 'temps lec', systime(1)-tempsun
352;
353   return,reform(z)
354
355end
356
Note: See TracBrowser for help on using the repository browser.