source: trunk/SRC/Obsolete/meshlec.pro @ 325

Last change on this file since 325 was 325, checked in by pinsard, 17 years ago

modification of some headers (+some corrections) to prepare usage of the new idldoc

  • Property svn:keywords set to Id
File size: 14.5 KB
Line 
1;+
2;
3; @file_comments
4; lecture du mask des sorties d'OPA. les sources se trouvent ds les
5; repertoires sur maia du type:
6; /nom_exp/RESTARTS
7;
8; @obsolete
9;
10; @examples
11; IDL> meshmask[,' nomfich']
12;
13; @param nomfich {in}{required}
14; string, c''est le nom du fichier a lire. Par defaut, c''est meshmask.
15;
16; @keyword GLAMBOUNDARY {in}
17; un vecteur de 2 elements specifaint le min et le
18; max qui doivent etre imposes en longitude (obligatoire si le
19; tableau depasse 360 degres)
20;
21; @keyword PASBLABLA {in}
22; pour supprimer les blablas
23;
24; @keyword DOUBLE {in}
25; pour forcer a lire les tableaux en double
26;        precision. ce Mot clef est maintenant active
27;        automatiquement.
28;
29; @keyword GETDIMENSIONS
30;
31; @uses
32; common.pro
33;
34; @restrictions
35;  La definition de ixminmesh,ixmaxmesh,iyminmesh,iymaxmesh
36;  ,izminmesh,izmaxmesh doit etre faite avant l''entree dans cette
37;  routine. pour attribuer automatiquement ces valeurs au maximum
38;  possible les mettre toutes a -1 et meshlec les calculera.
39;
40; @history
41; Sebastien Masson (smasson\@lodyc.jussieu.fr)
42;     Marina Levy : lecture en double precision (cas calcul sur shine)
43;
44; @version
45; $Id$
46;
47;-
48PRO meshlec, nomfich, PASBLABLA = pasblabla,DOUBLE = double, GLAMBOUNDARY = glamboundary, GETDIMENSIONS = GETDIMENSIONS
49;
50  compile_opt idl2, strictarrsubs, obsolete
51;
52@common
53   tempsun = systime(1)         ; pour key_performance
54;---------------------------------------------------------
55;---------------------------------------------------------
56   jpiglo      = 0L
57   jpjglo      = 0L
58   jpkglo      = 0L
59   tab    = 'aaaaa'
60;---------------------------------------------------------
61;  definition du domaine de la grille sur lequel sont
62;  effectuees les sorties.les indices des tableaux commencant
63;  a 1: cf le fichier wrivr2.F ds WKOPA sur le cray
64;----------------------------------------------------------
65; LECTURE DU MASK trouve ds les fichiers restart
66;----------------------------------------------------------
67;----------------------------------------------------------
68;
69;----------------------------------------------------------
70;       constitution de l'adresse s_fichier et
71;     ouverture du fichier a l'adresse s_fichier
72;-------------------------------------------------------
73
74   IF n_params() EQ 0 then nomfich = 'meshmask'
75   s_fichier = isafile(file = nomfich, iodir = iodir)
76   if not keyword_set(pasblabla) then print, ' '
77   if not keyword_set(pasblabla) then print,'adresse du fichier: ',s_fichier
78
79   openr, numlec, s_fichier, /get_lun, /f77_unformatted, /swap_if_little_endian
80   filepamameters = fstat(numlec)
81;---------------------------------------------------------
82;     lecture
83;---------------------------------------------------------
84   readu, numlec, jpiglo, jpjglo, jpkglo
85   if not keyword_set(pasblabla) then print, 'taille de la grille d''origine: ',jpiglo,', ',jpjglo,', ',jpkglo
86;
87   if keyword_set(getdimensions) then begin
88      free_lun,numlec
89      close, numlec
90      return
91   endif
92
93;---------------------------------------------------------
94; on determine si le fichier a ete ecrit en double precision on non...
95   sizenumber = 8l
96   sizefile8 = 4l+3l*4l+4l $
97    +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
98    +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
99    +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
100    +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
101    +4l*(4l+jpiglo*jpjglo*jpkglo*sizenumber+4l) $
102    +1l*(4l+jpiglo*jpjglo*sizenumber+4l) $
103    +4l*(4l+jpkglo*sizenumber+4l)
104   if filepamameters.size GE sizefile8 THEN double = 1
105;    sizenumber = 4l
106;    sizefile4 = 4l+3l*4l+4l $
107;     +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
108;     +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
109;     +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
110;     +4l*(4l+jpiglo*jpjglo*sizenumber+4l) $
111;     +4l*(4l+jpiglo*jpjglo*jpkglo*sizenumber+4l) $
112;     +1l*(4l+jpiglo*jpjglo*sizenumber+4l) $
113;     +4l*(4l+jpkglo*sizenumber+4l)
114; print, filepamameters.size,  sizefile4,  sizefile8
115;    case filepamameters.size of
116;       sizefile8:double = 1
117;       sizefile4:double = 0
118;       ELSE:BEGIN
119;          nothing = report('The OPA Mesh file as not the good size!')
120;          free_lun,numlec
121;          close, numlec
122;          return
123;       END
124;    endcase
125;
126;
127   if n_elements(ixminmesh) EQ 0 THEN ixminmesh = 0
128   if n_elements(ixmaxmesh) EQ 0 then ixmaxmesh = jpiglo-1
129   if ixminmesh EQ -1 THEN ixminmesh = 0
130   IF ixmaxmesh EQ -1 then ixmaxmesh = jpiglo-1
131   if n_elements(iyminmesh) EQ 0 THEN iyminmesh = 0
132   IF n_elements(iymaxmesh) EQ 0 then iymaxmesh = jpjglo-1
133   if iyminmesh EQ -1 THEN iyminmesh = 0
134   IF iymaxmesh EQ -1 then iymaxmesh = jpjglo-1
135   if n_elements(izminmesh) EQ 0 THEN izminmesh = 0
136   IF n_elements(izmaxmesh) EQ 0 then izmaxmesh = jpkglo-1
137   if izminmesh EQ -1 THEN izminmesh = 0
138   IF izmaxmesh EQ -1 then izmaxmesh = jpkglo-1
139;
140;
141
142   jpi    = long(ixmaxmesh-ixminmesh+1)
143   jpj    = long(iymaxmesh-iyminmesh+1)
144   jpk    = long(izmaxmesh-izminmesh+1)
145;-------------------------------------------------------
146; doit-on reellement lire la grille?
147;-------------------------------------------------------
148   meshparameters = {jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $
149                     , jpi:jpi, jpj:jpj, jpk:jpk $
150                     , ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $
151                     , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $
152                     , izminmesh:izminmesh, izmaxmesh:izmaxmesh $
153                     , key_shift:key_shift}
154;
155;-------------------------------------------------------
156   noticebase=xnotice('Lecture du fichier !C '+s_fichier+'!C ...')
157;-------------------------------------------------------
158
159   IF NOT keyword_set(double) THEN BEGIN
160      z3d    = fltarr(jpiglo, jpjglo, jpkglo)
161      z2d    = fltarr(jpiglo, jpjglo)
162      z1d    = fltarr(jpkglo)
163   ENDIF ELSE BEGIN
164      z3d    = dblarr(jpiglo, jpjglo, jpkglo)
165      z2d    = dblarr(jpiglo, jpjglo)
166      z1d    = dblarr(jpkglo)
167   ENDELSE
168
169   if not keyword_set(pasblabla) then print, ' '
170   readu, numlec, tab,z2d
171   GLAMT=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
172   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GLAMT[25,31]: ',GLAMT[25,31]
173   readu, numlec, tab,z2d
174   GLAMU=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
175   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GLAMU[25,31]: ',GLAMU[25,31]
176   readu, numlec, tab,z2d
177   GLAMV=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
178   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GLAMV[25,31]: ',GLAMV[25,31]
179   readu, numlec, tab,z2d
180   GLAMF=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
181   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GLAMF[25,31]: ',z2d[25,31]
182
183   if not keyword_set(pasblabla) then print, ' '
184   readu, numlec, tab, z2d
185   GPHIT=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
186   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GPHIT[25,31]: ',GPHIT[25,31]
187   readu, numlec, tab, z2d
188   GPHIU=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
189   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GPHIU[25,31]: ',GPHIU[25,31]
190   readu, numlec, tab, z2d
191   GPHIV=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
192   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GPHIV[25,31]: ',GPHIV[25,31]
193   readu, numlec, tab, z2d
194   GPHIF=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
195   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GPHIF[25,31]: ',z2d[25,31]
196
197   if not keyword_set(pasblabla) then print, ' '
198   readu, numlec, tab, z2d
199   E1T=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
200   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E1T[25,5]: ',z2d[25,5]
201   readu, numlec, tab, z2d
202   E1U=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
203   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E1U[25,5]: ',z2d[25,5]
204   readu, numlec, tab, z2d
205   E1V=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
206   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E1V[25,5]: ',z2d[25,5]
207   readu, numlec, tab, z2d
208   E1F=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
209   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E1F[25,5]: ',z2d[25,5]
210
211   if not keyword_set(pasblabla) then print, ' '
212   readu, numlec, tab, z2d
213   E2T=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
214   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E2T[25,5]: ',z2d[25,5]
215   readu, numlec, tab, z2d
216   E2U=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
217   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E2U[25,5]: ',z2d[25,5]
218   readu, numlec, tab, z2d
219   E2V=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
220   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E2V[25,5]: ',z2d[25,5]
221   readu, numlec, tab, z2d
222   E2F=float(z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh])
223   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E2F[25,5]: ',z2d[25,5]
224
225   if not keyword_set(pasblabla) then print, ' '
226   readu, numlec, tab, z3d
227   TMASK=byte(z3d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh,izminmesh:izmaxmesh])
228   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur TMASK[25,5,0]: ',TMASK[25,5,0]
229   readu, numlec, tab, z3d
230   UMASKred=byte(z3d[ixmaxmesh,iyminmesh:iymaxmesh,izminmesh:izmaxmesh])
231   umaskred = reform(umaskred)
232   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur UMASK[25,5,0]: ',z3d[25,5,0]
233   readu, numlec, tab, z3d
234   VMASKred=byte(z3d[ixminmesh:ixmaxmesh,iymaxmesh,izminmesh:izmaxmesh])
235   vmaskred = reform(vmaskred)
236   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur VMASK[25,5,0]: ',z3d[25,5,0]
237   readu, numlec, tab, z3d
238   fmaskredy=byte(z3d[ixmaxmesh,iyminmesh:iymaxmesh,izminmesh:izmaxmesh])
239   coast = where(fmaskredy NE 0 and fmaskredy NE 1)
240   IF coast[0] NE -1 THEN fmaskredy[coast] = 0b
241   fmaskredx=byte(z3d[ixminmesh:ixmaxmesh,iymaxmesh,izminmesh:izmaxmesh])
242   coast = where(fmaskredx NE 0 and fmaskredx NE 1)
243   IF coast[0] NE -1 THEN fmaskredx[coast] = 0b
244   fmaskredx = reform(fmaskredx)
245   fmaskredy = reform(fmaskredy)
246   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur FMASK[25,5,0]: ',z3d[25,5,0]
247;
248   if not keyword_set(pasblabla) then print, ' '
249   readu, numlec, tab, z2d
250;FF=z2d[ixminmesh:ixmaxmesh,iyminmesh:iymaxmesh]
251   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur FF[25,5]: ',z2d[25,5]
252   readu, numlec, tab, z1d
253   GDEPT=float(z1d[izminmesh:izmaxmesh])
254   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GDEPT[1]: ',GDEPT[1]
255   readu, numlec, tab, z1d
256   GDEPW=float(z1d[izminmesh:izmaxmesh])
257   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur GDEPW[1]: ',GDEPW[1]
258   readu, numlec, tab, z1d
259   E3T=float(z1d[izminmesh:izmaxmesh])
260   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E3T[3]: ',E3T[3]
261   readu, numlec, tab, z1d
262   E3W=float(z1d[izminmesh:izmaxmesh])
263   if not keyword_set(pasblabla) then print, 'tableau: ',tab,' exemple de valeur E3W[3]: ',E3W[3]
264;-------------------------------------------------------
265   free_lun,numlec
266   close, numlec
267;-------------------------------------------------------
268;-------------------------------------------------------
269; bornes de glam qui ne doivent pas depasser 360 degres..
270;-------------------------------------------------------
271;    minglam = min(glamt, max = maxglam)
272;    if maxglam-minglam GE 360 AND NOT keyword_set(glamboundary) then $
273;     nothing = execute('glamboundary = '+xquestion('What are the longitudes boundary?', '[-180,180]',/chkwidget))
274;-------------------------------------------------------
275   if keyword_set(glamboundary) then begin
276      if glamboundary[0] NE glamboundary[1] then begin
277         glamt = glamt MOD 360
278         smaller = where(glamt LT glamboundary[0])
279         if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360
280         bigger = where(glamt GE glamboundary[1])
281         if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360
282         glamu = glamu MOD 360
283         smaller = where(glamu LT glamboundary[0])
284         if smaller[0] NE -1 then glamu[smaller] = glamu[smaller]+360
285         bigger = where(glamu GE glamboundary[1])
286         if bigger[0] NE -1 then glamu[bigger] = glamu[bigger]-360
287         glamv = glamv MOD 360
288         smaller = where(glamv LT glamboundary[0])
289         if smaller[0] NE -1 then glamv[smaller] = glamv[smaller]+360
290         bigger = where(glamv GE glamboundary[1])
291         if bigger[0] NE -1 then glamv[bigger] = glamv[bigger]-360
292         glamf = glamf MOD 360
293         smaller = where(glamf LT glamboundary[0])
294         if smaller[0] NE -1 then glamf[smaller] = glamf[smaller]+360
295         bigger = where(glamf GE glamboundary[1])
296         if bigger[0] NE -1 then glamf[bigger] = glamf[bigger]-360
297      endif
298   endif
299;-------------------------------------------------------
300; shift en x
301;-------------------------------------------------------
302   if keyword_set(key_shift) AND jpi NE 1 then begin
303      glamt = shift(glamt,key_shift , 0)
304      gphit = shift(gphit,key_shift , 0)
305      e1t = shift(e1t, key_shift, 0)
306      e2t = shift(e2t,key_shift , 0)
307      glamu = shift(glamu, key_shift, 0)
308      gphiu = shift(gphiu, key_shift, 0)
309      e1u = shift(e1u,key_shift , 0)
310      e2u = shift(e2u, key_shift, 0)
311      glamv = shift(glamv, key_shift, 0)
312      gphiv = shift(gphiv, key_shift, 0)
313      e1v = shift(e1v,key_shift , 0)
314      e2v = shift(e2v, key_shift, 0)
315      glamf = shift(glamf, key_shift, 0)
316      gphif = shift(gphif, key_shift, 0)
317      e1f = shift(e1f, key_shift , 0)
318      e2f = shift(e2f, key_shift, 0)
319      if jpk EQ 1 then begin
320         tmask = shift(tmask, key_shift, 0)
321         vmaskred = shift(vmaskred, key_shift)
322         fmaskredx = shift(fmaskredx, key_shift)
323      ENDIF ELSE BEGIN
324         tmask = shift(tmask, key_shift, 0,0)
325         vmaskred = shift(vmaskred, key_shift, 0)
326         fmaskredx = shift(fmaskredx, key_shift, 0)
327      ENDELSE
328   endif
329;-------------------------------------------------------
330   key_yreverse = 0
331   key_zreverse = 0
332   key_partialstep = 0
333   key_stride = [1, 1, 1]
334   key_gridtype = 'c'
335;
336   if not keyword_set(pasblabla) then print,'lecture '+nomfich+' finie'
337   widget_control, noticebase, bad_id = toto, /destroy
338   if keyword_set(key_performance) THEN print, 'temps meshlec', systime(1)-tempsun
339
340   return
341end
342
343
344
345
Note: See TracBrowser for help on using the repository browser.