source: trunk/condmag_from_orca.pro @ 21

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

some improvements on validation of step 1 and begining of parametrisation to orcares
  • Property svn:keywords set to Id
File size: 13.0 KB
Line 
1;+
2;
3; @file_comments
4; interpolate cond_sed_<orcares>.nc and Br_<orcares>.nc
5; and produce condmag_from_<orcares>.nc file
6;
7; condmag_from_<orcares>.nc might be compare to condmag.nc
8;
9; NB : cond_oc and cond_tot won't be reproduced
10; ++ should I write missing values ?
11;
12; inverse of <progeomag>condmag_on_orca</progeomag>
13;
14; @categories
15; interpolation, orca grid
16;
17; @param orcares {in}{required}{type=string}
18; code of ORCA grid to use for output results
19; must be 'ORCA2' or 'ORCA025'
20;
21; @param method {in}{required}{type=string}
22; code for interpolation method
23; must be 'bilinear'
24; parameter is therefore here to be visible. see restrictions below.
25;
26; @param gridtype {in}{required}{type=string}
27; to specify on which grid we do the interpolation T, U, V
28; must belong to T, U or V
29; ++ je ne sais pas comment ce paramètre doit intervenir
30;
31; @keyword DRAKKAR_EXP {type=string}
32; code for Drakkar experiment
33; only used when orcares = ORCA025
34; must be G42 ++ G70
35; ++ je ne sais pas comment ce paramètre doit intervenir a priori
36; la partie mesh est la même pour toutes les expériences mais attention
37; pas forcément le mask.
38;
39; @keyword PERF
40; if set activate performance with <prodidl>profiler</proidl>
41; PERF vs key_performance dans les commons
42;
43; @restrictions
44;  - Requires <a href="http://forge.ipsl.jussieu.fr/saxo/">SAXO</a>
45;  - won't be relevant if 'imoms3' was used as interpolation method in
46; <progeomag>condmag_on_orca</progeomag>
47;
48; @todo
49; add version variable inside code
50; provide tools to plot output files
51; produce a NetCDF GDT or CF compliant
52; introduce <progeomag>../<progeomag> and <pro>...</pro> in postprocessing
53; of idldoc outputs
54;
55; @pre
56; see geomag_profile.sh
57;
58; be sure to have cond_sed_*.nc and Br_*.nc in the directory defined
59; in ${GEOMAG_ID}/
60; be sure to have condmag.nc in the directory defined
61; in ${GEOMAG_ID}/
62; be sure to not have condmag_from_*.nc in the directory defined in
63; ${GEOMAG_OD}/
64; be sure to have meshmask of ORCA grid you choose in the directory defined in
65; ${GEOMAG_ID}/
66; for ORCA2 filename is meshmask_bab.nc
67; ++ au pif entre mesh_hgr, mesh_zgr et mask
68; for ORCA025 filename is ORCA025-G42_mesh_hgr.nc
69;
70; @post
71; see geomag_profile.sh
72; condmag_from_orcares.nc is now present in ${GEOMAG_OD}/
73;
74; @examples
75; IDL> .run condmag_from_orca
76; IDL> condmag_from_orca, 'ORCA2', 'bilinear','T'
77;
78; @history
79; fplod 2007-03-21T13:23:55Z aedon.locean-ipsl.upmc.fr (Darwin)
80; create from condmag_on_orca to validate step1 of GEOMAG
81;
82; @version
83; $Id$
84;
85;-
86PRO condmag_from_orca, orcares, method, gridtype, DRAKKAR_EXP = drakkar_exp, PERF = perf
87;
88  compile_opt idl2, strictarrsubs
89;
90IF keyword_set(perf) EQ 1 THEN BEGIN
91    msg = 'iii : start profiler'
92    ras = report(msg)
93    ;PROFILER, /SYSTEM & PROFILER
94    PROFILER, /SYSTEM & PROFILER
95; ++ ne tracera que les modules compilés au moment de l'appel à cette commande
96; donc ++ regarder si on a bien tout et si on doit ajouter des  PROFILER,toto dans
97; tous les modules
98ENDIF
99;
100;----
101; check input parameters
102;----
103;
104; check orcares definition
105;
106  CASE orcares OF
107     'ORCA2': BEGIN
108                 msg = 'iii : valid orcares parameter = ' + orcares
109                 ras = report(msg)
110                 filename_oce = 'meshmask_bab.nc'
111                 IF keyword_set(DRAKKAR_EXP) THEN BEGIN
112                    msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp
113                    ras = report(msg)
114                 END
115              END
116     'ORCA025': BEGIN
117                 msg = 'iii : valid orcares parameter = ' + orcares
118                 ras = report(msg)
119                 IF keyword_set(DRAKKAR_EXP) THEN BEGIN
120                    msg = 'iii : DRAKKAR_EXP keyword set'
121                    ras = report(msg)
122                    msg = 'iii : DRAKKAR_EXP = ' + drakkar_exp
123                    ras = report(msg)
124                    CASE drakkar_exp OF
125                       'G42' : BEGIN
126                                  msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp
127                                  ras = report(msg)
128                               END
129                       'G70' : BEGIN
130                                  msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp
131                                  ras = report(msg)
132                               END
133                       ELSE : BEGIN
134                                  msg = 'eee : invalid DRAKKAR_EXP keyword = ' + drakkar_exp
135                                  ras = report(msg)
136                                  RETURN
137                              END
138                    ENDCASE
139                    filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc'
140                 END
141              END
142      ELSE  : BEGIN
143                 msg = 'eee : invalid orcares parameter = ' + orcares
144                 ras = report(msg)
145                 msg = 'eee : orcares must be ORCA2 or ORCA025'
146                 ras = report(msg)
147                 RETURN
148              END
149  ENDCASE
150;
151; check method definition
152  CASE method OF
153     'bilinear': BEGIN
154                    msg = 'iii : valid method parameter = ' + method
155                    ras = report(msg)
156                 END
157     ELSE      : BEGIN
158                    msg = 'eee : invalid method parameter = ' + method
159                    ras = report(msg)
160                    msg = 'eee : method must be bilinear'
161                    ras = report(msg)
162                    RETURN
163                 END
164  ENDCASE
165;
166; check gridtype definition
167  gridtype = strupcase(gridtype)
168  CASE gridtype OF
169     'T' : BEGIN
170              msg = 'iii : valid gridtype parameter = ' + gridtype
171              ras = report(msg)
172           END
173     'U' : BEGIN
174              msg = 'iii : valid gridtype parameter = ' + gridtype
175              ras = report(msg)
176           END
177     'V' : BEGIN
178              msg = 'iii : valid gridtype parameter = ' + gridtype
179              ras = report(msg)
180           END
181     ELSE: BEGIN
182              msg = 'eee : invalid gridtype parameter = ' + gridtype
183              ras = report(msg)
184              msg = 'eee : gridtype must be T, U or V'
185              ras = report(msg)
186              RETURN
187           END
188  ENDCASE
189;
190; check for input files
191;
192; test if ${GEOMAG_ID} defined
193  geomag_id_env=GETENV('GEOMAG_ID')
194  CASE geomag_id_env OF
195     ''  :  BEGIN
196              msg = 'eee : ${GEOMAG_ID} is not defined'
197              ras = report(msg)
198              RETURN
199            END
200     ELSE: BEGIN
201             msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env
202             ras = report(msg)
203           END
204  ENDCASE
205;
206  iodirin = isadirectory(geomag_id_env)
207;
208; existence and protection of ${GEOMAG_ID}
209  IF (FILE_TEST(iodirin, /DIRECTORY,/EXECUTABLE , /READ) EQ 0) THEN BEGIN
210     msg = 'eee : the directory' + iodirin  + ' is not accessible.'
211     ras = report(msg)
212     RETURN
213  ENDIF
214;
215; build input filenames
216  filename_cond_sed = 'cond_sed' + '_' + orcares +'.nc'
217  fullfilename_cond_sed = iodirin + filename_cond_sed
218;
219; check if this file exists
220  fullfilename_cond_sed = isafile(iodirin + filename_cond_sed, NEW=0,/MUST_EXIST)
221  IF fullfilename_cond_sed[0] EQ '' THEN BEGIN
222     msg = 'eee : the file ' + fullfilename_cond_sed + ' was not found.'
223     ras = report(msg)
224     RETURN
225  ENDIF
226;
227; protection
228  IF (FILE_TEST(fullfilename_cond_sed[0], /READ) EQ 0) THEN BEGIN
229     msg = 'eee : the file ' + fullfilename_cond_sed[0] + ' is not readable.'
230     ras = report(msg)
231     RETURN
232  ENDIF
233;
234  filename_Br = 'Br' + '_' + orcares +'.nc'
235  fullfilename_Br = iodirin + filename_Br
236;
237; check if this file exists
238  fullfilename_Br = isafile(iodirin + filename_Br, NEW=0,/MUST_EXIST)
239  IF fullfilename_Br[0] EQ '' THEN BEGIN
240     msg = 'eee : the file ' + fullfilename_Br + ' was not found.'
241     ras = report(msg)
242     RETURN
243  ENDIF
244;
245; protection
246  IF (FILE_TEST(fullfilename_Br[0], /READ) EQ 0) THEN BEGIN
247     msg = 'eee : the file ' + fullfilename_Br[0] + ' is not readable.'
248     ras = report(msg)
249     RETURN
250  ENDIF
251;
252; mesh mask
253; check if this file exists
254  fullfilename_oce = isafile(iodirin + filename_oce, NEW=0, /MUST_EXIST, $
255                     RECURSIVE=0)
256  IF fullfilename_oce[0] EQ '' THEN BEGIN
257     msg = 'eee : the file ' + fullfilename_oce + ' was not found.'
258     ras = report(msg)
259     RETURN
260  ENDIF
261;
262; protection
263  IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN
264     msg = 'eee : the file ' + fullfilename_oce[0] + ' is not readable.'
265     ras = report(msg)
266     RETURN
267  ENDIF
268;
269  filename_condmag = 'condmag.nc'
270  fullfilename_condmag = iodirin + filename_condmag
271;
272; check if this file exists
273  fullfilename_condmag = isafile(iodirin + filename_condmag, NEW=0,/MUST_EXIST)
274  IF fullfilename_condmag[0] EQ '' THEN BEGIN
275     msg = 'eee : the file ' + fullfilename_condmag + ' was not found.'
276     ras = report(msg)
277     RETURN
278  ENDIF
279;
280; protection
281  IF (FILE_TEST(fullfilename_condmag[0], /READ) EQ 0) THEN BEGIN
282     msg = 'eee : the file ' + fullfilename_condmag[0] + ' is not readable.'
283     ras = report(msg)
284     RETURN
285  ENDIF
286;
287; test if ${GEOMAG_OD} defined
288  geomag_od_env=GETENV('GEOMAG_OD')
289  CASE geomag_od_env OF
290     '' : BEGIN
291             msg = 'eee : ${GEOMAG_OD} is not defined'
292             ras = report(msg)
293             RETURN
294          END
295     ELSE: BEGIN
296             msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env
297             ras = report(msg)
298           END
299  ENDCASE
300;
301; check if output data will be possible
302  iodirout = isadirectory(geomag_od_env)
303;
304; existence and protection
305  IF (FILE_TEST(iodirout, /DIRECTORY,/WRITE) EQ 0) THEN BEGIN
306     msg = 'eee : the directory' + iodirout  + ' was not found.'
307     ras = report(msg)
308     RETURN
309  ENDIF
310;
311; build output filename
312  filename = 'condmag_from_' +  orcares +'.nc'
313  fullfilename = iodirout + filename
314;
315; in order to avoid unexpected overwritten
316  IF (FILE_TEST(fullfilename) EQ 1) THEN BEGIN
317     msg = 'eee : the file ' + fullfilename  + ' already exists.'
318     ras = report(msg)
319     RETURN
320  ENDIF
321;
322;----
323; Oceanic grid parameters ie input grid + mask
324;----
325;
326CASE orcares OF
327  'ORCA2': BEGIN
328              @initorca2_bab
329           END
330  ENDCASE
331;
332;----
333; read data to interpolate on regular grid
334;----
335;
336a=read_ncdf('cond_sed', 0, /timestep, $
337file = fullfilename_cond_sed[0], /nostruct)
338msg = 'iii : ' + fullfilename_cond_sed[0] + ' opened for read'
339ras = report(msg)
340; mask and limite cond_sed on values below 5000 ++
341cond_sedin = extrapsmooth(a, tmask[*, *, 0]*(a le 5000), /x_periodic)
342;
343b = read_ncdf('Br', 0, /timestep, $
344file = fullfilename_Br[0], /nostruct)
345msg = 'iii : ' + fullfilename_Br[0] + ' opened for read'
346ras = report(msg)
347; mask and limit Br on values below/above ++
348Brin = extrapsmooth(b, tmask[*, *, 0]*(b), /x_periodic)
349;
350glamin = glamt
351gphiin = gphit
352maskin = tmask[*, *, 0]
353;
354;---
355; conductivity and magnetic field grid parameters ie output grid
356;----
357
358initncdf, fullfilename_condmag[0], xaxisname = 'lo', yaxisname = 'la'
359msg = 'iii : ' + fullfilename_condmag[0] + ' opened for read'
360ras = report(msg)
361;
362glamout = glamt
363gphiout = gphit
364;
365
366dimidxout = jpi
367dimidyout = jpj
368;
369;---------------
370; Interpolate
371;---------------
372;++ cond_sedout = fromirr('bilinear', cond_sedin, glamin, gphiin, maskin,  glamout, gphiout, -1,  WEIG = weights, ADDR = addresses)
373cond_sedout=dindgen(dimidxout,dimidyout)
374;Brout=fromirr('bilinear', Brin, weights, addresses) ; ++ pb doc % Variable is undefined: LONOUT
375;++Brout = fromirr('bilinear', Brin, glamin, gphiin, maskin,  glamout, gphiout, -1,  WEIG = weights, ADDR = addresses)
376;
377Brout=dindgen(dimidxout,dimidyout)
378
379;
380;---------------
381; Produce outputs
382;---------------
383;
384  netcdf_id = NCDF_CREATE(fullfilename, /clobber)
385  NCDF_CONTROL, netcdf_id, /NOFILL
386;
387; dimension
388  dimidxout = NCDF_DIMDEF(netcdf_id, 'la' ,  jpi)
389  dimidyout = NCDF_DIMDEF(netcdf_id, 'lo' ,  jpj)
390;
391; global attributes
392  NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite
393;++  NCDF_ATTPUT, netcdf_id, 'file_name'  , filename, /GLOBAL
394;++  NCDF_ATTPUT, netcdf_id, 'Title'      , title, /GLOBAL
395;
396; declaration of variables
397  varid = lonarr(4)
398;
399  varid[0] = NCDF_VARDEF(netcdf_id, 'la'  , [dimidxout], /FLOAT)
400  varid[1] = NCDF_VARDEF(netcdf_id, 'lo'  , [dimidyout], /FLOAT)
401  varid[2] = NCDF_VARDEF(netcdf_id, 'cond_sed'  , [dimidxout, dimidyout], /FLOAT)
402  NCDF_ATTPUT, netcdf_id, varid[2], 'units'    , 'siemens'
403  varid[3] = NCDF_VARDEF(netcdf_id, 'Br'  , [dimidxout, dimidyout], /FLOAT)
404  NCDF_ATTPUT, netcdf_id, varid[3], 'units'    , 'siemens'
405;
406  NCDF_CONTROL, netcdf_id, /ENDEF
407;
408  NCDF_VARPUT, netcdf_id, varid[0], dindgen(dimidxout,dimidyout)
409  NCDF_VARPUT, netcdf_id, varid[1], dindgen(dimidyout,dimidyout)
410;
411; write the data
412;
413  NCDF_VARPUT, netcdf_id, varid[2], cond_sedout
414;++  NCDF_VARPUT, netcdf_id, 'Br', Brout
415;---------------------------
416; close the netcdf files
417  NCDF_CLOSE, netcdf_id
418;
419  msg = 'iii : ' + fullfilename + ' created'
420  ras = report(msg)
421;
422IF keyword_set(perf) EQ 1 THEN BEGIN
423  msg = 'iii : report profiler results'
424  ras = report(msg)
425; ++ tri par ordre alpha , par ordre de pourcentage croissant ou décroissant
426; +++ d'utilisation
427  profiler,/REPORT
428;  shut down all profiling (according to
429; http://www.dfanning.com/code_tips/whyslow.html)
430  profiler, /CLEAR, /SYSTEM
431  profiler, /CLEAR, /RESET
432ENDIF
433;
434END
Note: See TracBrowser for help on using the repository browser.