source: trunk/GRILLE/ncdf_meshlec.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: 11.9 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:ncdf_meshlec
6;
7; PURPOSE:lit le fichier meshmask au format netcdf cree par OPA ds sa
8; version tout NETCDF!
9;
10; CATEGORY:lecture de grille
11;
12; CALLING SEQUENCE:ncdf_meshlec [,' nomfich']
13;
14; INPUTS:
15;
16;       nomfich: un string definissant le nom du fichier a lire. Si
17; ce string ne contient pas l''adresse entiere du fichier nomfich est
18; complete par iodir. Par defaut nomfich=iodir+'meshmask.nc'
19;
20; KEYWORD PARAMETERS:
21;
22;    GLAMBOUNDARY: un vecteur de 2 elements specifaint le min et le
23;    max qui doivent etre imposes en longitude (obligatoire si le
24;    tableau depasse 360 degres)
25;
26;    /CHECKDAT: cherche si il existe un ficher du meme nom que le
27;    fichier meshmask mais terminant par .dat au lieu de .dat.
28;             Si ce fichier n''existe pas, ncdf_meshlec lit le fichier
29;         NetCdf et sauve tous les tableaux lus ds le fichier
30;         ...dat. Rq: Ce fichier est bcp plus petit que le fichier
31;         NetCdf d''origine car par ex on ne sauve pas les tableaux
32;         umask, vmask, fmask mais simplement certains de leurs bords
33;         (cf. umask.pro vmask.pro et fmask.pro) et ces masks sont
34;         sauves en bytes et toues les autres tableaux en
35;         float. Cependant ce fichier .dat est cree en fonction des
36;         parametres rentres ds le init... (ixminmesh.....) et Je ne
37;         sais rien sur la portabilite du fichier .dat.
38;            Si ce fichier .dat existe on le lit avec restore.
39;
40; OUTPUTS:mise a jours des variables de common glam, common
41; gphi,common e12,common mask,common tab (cf. common.pro)
42;
43; COMMON BLOCKS:common.pro
44;
45; SIDE EFFECTS:prend en compte (si il sont definits) les
46; ixminmesh,...et definit jpiglo,jpi,....
47;
48; RESTRICTIONS:
49;
50;  La definition de ixminmesh,ixmaxmesh,iyminmesh,iymaxmesh
51;  ,izminmesh,izmaxmesh doit etre faite avant l''entree dans cette
52;  routine. pour attribuer automatiquement ces valeurs au maximum
53;  possible les mettre toutes a -1 et ncdf_meshlec les calculera.
54;
55; EXAMPLE:
56;
57; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
58;                      12/1999
59;-
60;------------------------------------------------------------
61;------------------------------------------------------------
62;------------------------------------------------------------
63PRO ncdf_meshlec, nomfich, GLAMBOUNDARY = glamboundary, GETDIMENSIONS = getdimensions, CHECKDAT = checkdat
64@common
65   tempsun = systime(1)         ; pour key_performance
66;----------------------------------------------------------
67;       constitution de l''adresse s_fichier et
68;     ouverture du fichier a l''adresse s_fichier
69;-------------------------------------------------------
70; def de nomfich par defaut
71   IF n_params() EQ 0 then nomfich = 'meshmask.nc'
72;------------------
73   if keyword_set(CHECKDAT) then begin
74      nomfichdat = strmid(nomfich, 0, strlen(nomfich)-2)+'dat'
75      thisOS = strupcase(strmid(!version.os_family, 0, 3))
76      CASE thisOS OF
77         'MAC':sep = ':'
78         'WIN':sep = '\'
79         ELSE:sep = '/'
80      ENDCASE
81      if strpos(nomfichdat, sep) EQ -1 then begin
82; si iodir ne finit pas par sep on le rajoute
83         if rstrpos(iodir, sep) NE strlen(iodir)-1 then iodir = iodir+sep
84         nomfichdat = iodir+nomfichdat
85      endif
86      test = findfile(nomfichdat)  ; le nom cherche correspond bien a un fichier?
87      if test[0] NE '' then begin
88         noticebase=xnotice('Lecture du fichier !C '+nomfichdat+'!C ...')
89         restore, nomfichdat
90         widget_control, noticebase, bad_id = nothing, /destroy
91         if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun
92         return
93      ENDIF
94   endif
95;------------------
96   s_fichier = isafile(file = nomfich, iodir = iodir)
97;
98   noticebase=xnotice('Lecture du fichier !C '+s_fichier+'!C ...')
99;
100   cdfid=ncdf_open(s_fichier)
101   contient=ncdf_inquire(cdfid)
102
103;------------------------------------------------------------
104; differentes dimensions
105;------------------------------------------------------------
106   ncdf_diminq,cdfid,'x',name,jpiglo
107   ncdf_diminq,cdfid,'y',name,jpjglo
108   ncdf_diminq,cdfid,'z',name,jpkglo
109;
110   if keyword_set(getdimensions) then begin
111      ncdf_close,  cdfid
112      return
113   endif
114;-------------------------------------------------------
115;-------------------------------------------------------
116   if n_elements(ixminmesh) EQ 0 THEN ixminmesh = 0
117   if n_elements(ixmaxmesh) EQ 0 then ixmaxmesh = jpiglo-1
118   if ixminmesh EQ -1 THEN ixminmesh = 0
119   IF ixmaxmesh EQ -1 then ixmaxmesh = jpiglo-1
120   if n_elements(iyminmesh) EQ 0 THEN iyminmesh = 0
121   IF n_elements(iymaxmesh) EQ 0 then iymaxmesh = jpjglo-1
122   if iyminmesh EQ -1 THEN iyminmesh = 0
123   IF iymaxmesh EQ -1 then iymaxmesh = jpjglo-1
124   if n_elements(izminmesh) EQ 0 THEN izminmesh = 0
125   IF n_elements(izmaxmesh) EQ 0 then izmaxmesh = jpkglo-1
126   if izminmesh EQ -1 THEN izminmesh = 0
127   IF izmaxmesh EQ -1 then izmaxmesh = jpkglo-1
128;
129   ixmindtasauve = testvar(var = ixmindta)
130   iymindtasauve = testvar(var = iymindta)
131   izmindtasauve = testvar(var = izmindta)
132;
133   ixmindta = 0l
134   iymindta = 0l
135   izmindta = 0l
136;
137   jpi    = long(ixmaxmesh-ixminmesh+1)
138   jpj    = long(iymaxmesh-iyminmesh+1)
139   jpk    = long(izmaxmesh-izminmesh+1)
140;
141   if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1]
142   key_stride = 1l > long(key_stride)
143;
144   jpitotal = jpi
145   key_shift = long(testvar(var = key_shift))
146   key = long(key_shift MOD jpitotal)
147   if key LT 0 then key = key+jpitotal
148   ixmin = ixminmesh-ixmindta
149;
150   jpi = jpi/key_stride[0]
151   jpj = jpj/key_stride[1]
152   jpk = jpk/key_stride[2]
153;
154   jpt = 1
155   premiertps = 0
156;
157   premierx = 0
158   dernierx = jpi-1
159   premiery = 0
160   derniery = jpj-1
161   premierz = 0
162   dernierz = jpk-1
163   nx = jpi
164   ny = jpj
165   nz = 1
166   izminmeshsauve = izminmesh
167   izminmesh = 0
168;
169;-------------------------------------------------------
170; tableaux 2d:
171;-------------------------------------------------------
172   if total(key_stride) EQ 3 then $
173    namevar = ['glamt','glamu','glamv','glamf','gphit','gphiu','gphiv','gphif','e1t','e1u','e1v','e1f','e2t','e2u','e2v','e2f'] $
174   ELSE namevar = ['glamt','glamu','glamv','gphit','gphiu','gphiv']
175;
176   for i = 0, n_elements(namevar)-1 do begin
177      varcontient=ncdf_varinq(cdfid, namevar[i])
178      name = varcontient.name
179@read_ncdf_varget
180      commande = namevar[i]+'=float(res)'
181      rien = execute(commande)
182   ENDFOR
183;
184   if total(key_stride) NE 3 then BEGIN
185      glamf = glamt
186      glamf = [glamf, reform(2.*reform(glamt[jpi-1, *])-reform(glamt[jpi-2, *]), 1, jpj)]
187      glamf = [ [glamf], [reform(2*glamf[*, jpj-1]-glamf[*, jpj-2], jpi+1, 1)] ]
188      glamf = glamf+shift(glamf, -1, -1)
189      glamf = .5*glamf[0:jpi-1, 0:jpj-1]
190;
191      gphif = gphit
192      gphif = [gphif, reform(2.*reform(gphit[jpi-1, *])-reform(gphit[jpi-2, *]), 1, jpj)]
193      gphif = [ [gphif], [reform(2*gphif[*, jpj-1]-gphif[*, jpj-2], jpi+1, 1)] ]
194      gphif = gphif+shift(gphif, -1, -1)
195      gphif = .5*gphif[0:jpi-1, 0:jpj-1]
196;
197      e1t = replicate(1, jpi, jpj)
198      e2t = e1t
199      e1u = e1t
200      e2u = e1t
201      e1v = e1t
202      e2v = e1t
203      e1f = e1t
204      e2f = e1t
205   ENDIF
206;-------------------------------------------------------
207; tableaux 3d:
208;-------------------------------------------------------
209   nz = jpk
210   izminmesh = izminmeshsauve
211;
212   varcontient=ncdf_varinq(cdfid,'tmask')
213   name = varcontient.name
214@read_ncdf_varget
215   tmask = byte(res)
216;
217   varcontient=ncdf_varinq(cdfid,'umask')
218   name = varcontient.name
219   nx = 1
220   premierx = jpi-1
221@read_ncdf_varget
222   umaskred = byte(res)
223   umaskred = reform(umaskred, /over)
224;
225   varcontient=ncdf_varinq(cdfid,'vmask')
226   name = varcontient.name
227   nx = jpi
228   premierx = 0
229   ny = 1
230   premiery = jpj-1
231@read_ncdf_varget
232   vmaskred = byte(res)
233   vmaskred = reform(vmaskred, /over)
234;
235   varcontient=ncdf_varinq(cdfid,'fmask')
236   name = varcontient.name
237   nx = 1
238   premierx = jpi-1
239   ny = jpj
240   premiery = 0
241@read_ncdf_varget
242   fmaskredy = byte(res)
243   coast = where(fmaskredy NE 0 and fmaskredy NE 1)
244   IF coast[0] NE -1 THEN fmaskredy[coast] = 0b
245   fmaskredy= reform(fmaskredy, /over)
246   nx = jpi
247   premierx = 0
248   ny = 1
249   premiery = jpj-1
250@read_ncdf_varget
251   fmaskredx = byte(res)
252   coast = where(fmaskredx NE 0 and fmaskredx NE 1)
253   IF coast[0] NE -1 THEN fmaskredx[coast] = 0b
254   fmaskredx = reform(fmaskredx, /over)
255;-------------------------------------------------------
256; tableaux 1d:
257;-------------------------------------------------------
258   namevar = ['e3t','e3w','gdept','gdepw']
259   for i = 0, n_elements(namevar)-1 do begin
260      varcontient=ncdf_varinq(cdfid, namevar[i])
261      if n_elements(varcontient.dim) EQ 4 THEN BEGIN
262         commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $
263          +',offset = [0,0,izminmesh,0], count = [1,1,jpk,1]'
264         if key_stride[2] NE 1 then commande = commande+', stride=[1,1,key_stride[2],1]'
265      ENDIF ELSE BEGIN
266         commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $
267          +',offset = [izminmesh], count = [jpk]'
268         if key_stride[2] NE 1 then commande = commande+', stride=key_stride[2]'
269      ENDELSE
270      rien = execute(commande)
271      commande = namevar[i]+'=float('+namevar[i]+')'
272      rien = execute(commande)
273      commande = namevar[i]+' = reform('+namevar[i]+', /over)'
274      rien = execute(commande)
275   ENDFOR
276;-------------------------------------------------------
277   ncdf_close,  cdfid
278;-------------------------------------------------------
279; bornes de glam qui ne doivent pas depasser 360 degres...
280;-------------------------------------------------------
281;    minglam = min(glamt, max = maxglam)
282;    if maxglam-minglam GE 360 AND NOT keyword_set(glamboundary) then $
283;     nothing = execute('glamboundary = '+xquestion('What are the longitudes boundary?', '[-180,180]',/chkwidget))
284;-------------------------------------------------------
285   if keyword_set(glamboundary) then BEGIN
286      if glamboundary[0] NE glamboundary[1] then BEGIN
287         glamt = glamt MOD 360
288         smaller = where(glamt LT glamboundary[0])
289         if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360
290         bigger = where(glamt GE glamboundary[1])
291         if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360
292         glamu = glamu MOD 360
293         smaller = where(glamu LT glamboundary[0])
294         if smaller[0] NE -1 then glamu[smaller] = glamu[smaller]+360
295         bigger = where(glamu GE glamboundary[1])
296         if bigger[0] NE -1 then glamu[bigger] = glamu[bigger]-360
297         glamv = glamv MOD 360
298         smaller = where(glamv LT glamboundary[0])
299         if smaller[0] NE -1 then glamv[smaller] = glamv[smaller]+360
300         bigger = where(glamv GE glamboundary[1])
301         if bigger[0] NE -1 then glamv[bigger] = glamv[bigger]-360
302         glamf = glamf MOD 360
303         smaller = where(glamf LT glamboundary[0])
304         if smaller[0] NE -1 then glamf[smaller] = glamf[smaller]+360
305         bigger = where(glamf GE glamboundary[1])
306         if bigger[0] NE -1 then glamf[bigger] = glamf[bigger]-360
307      endif
308   endif
309;-------------------------------------------------------
310   ixmindta = ixmindtasauve
311   iymindta = iymindtasauve
312   izmindta = izmindtasauve
313;-------------------------------------------------------
314   widget_control, noticebase, bad_id = nothing, /destroy
315;
316   if keyword_set(checkdat) then BEGIN
317   noticebase=xnotice('Ecriture du fichier !C '+nomfichdat+'!C ...')
318;
319      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
320;
321      widget_control, noticebase, bad_id = nothing, /destroy
322   endif
323;
324   if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun
325
326;-------------------------------------------------------
327   return
328end
Note: See TracBrowser for help on using the repository browser.