source: trunk/condmag_from_orca.pro @ 29

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

replace call of initorca2_bab by call of initocemesh
  • Property svn:keywords set to Id
File size: 13.4 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; reee522 2007-06-06T15:14:24Z rhodes (IRIX64)
80; replace call of initorca2_bab by initocemesh
81; reee522 2007-06-06T14:57:39Z rhodes (IRIX64)
82; correction for DRAKKAR_EXP usage
83; fplod 2007-03-21T13:23:55Z aedon.locean-ipsl.upmc.fr (Darwin)
84; create from condmag_on_orca to validate step1 of GEOMAG
85;
86; @version
87; $Id$
88;
89;-
90PRO condmag_from_orca, orcares, method, gridtype, DRAKKAR_EXP = drakkar_exp, PERF = perf
91;
92  compile_opt idl2, strictarrsubs
93;
94IF keyword_set(perf) EQ 1 THEN BEGIN
95    msg = 'iii : start profiler'
96    ras = report(msg)
97    ;PROFILER, /SYSTEM & PROFILER
98    PROFILER, /SYSTEM & PROFILER
99; ++ ne tracera que les modules compilés au moment de l'appel à cette commande
100; donc ++ regarder si on a bien tout et si on doit ajouter des  PROFILER,toto dans
101; tous les modules
102ENDIF
103;
104;----
105; check input parameters
106;----
107;
108; check orcares definition
109;
110  CASE orcares OF
111     'ORCA2': BEGIN
112                 msg = 'iii : valid orcares parameter = ' + orcares
113                 ras = report(msg)
114                 filename_oce = 'meshmask_bab.nc'
115                 IF keyword_set(DRAKKAR_EXP) THEN BEGIN
116                    msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp
117                    ras = report(msg)
118                 END
119              END
120     'ORCA025': BEGIN
121                 msg = 'iii : valid orcares parameter = ' + orcares
122                 ras = report(msg)
123                 IF keyword_set(DRAKKAR_EXP) THEN BEGIN
124                    msg = 'iii : DRAKKAR_EXP keyword set'
125                    ras = report(msg)
126                    msg = 'iii : DRAKKAR_EXP = ' + drakkar_exp
127                    ras = report(msg)
128                    CASE drakkar_exp OF
129                       'G42' : BEGIN
130                                  msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp
131                                  ras = report(msg)
132                               END
133                       'G70' : BEGIN
134                                  msg = 'iii : valid DRAKKAR_EXP keyword = ' + drakkar_exp
135                                  ras = report(msg)
136                               END
137                       ELSE : BEGIN
138                                  msg = 'eee : invalid DRAKKAR_EXP keyword = ' + drakkar_exp
139                                  ras = report(msg)
140                                  RETURN
141                              END
142                    ENDCASE
143                    filename_oce = orcares + '-' + drakkar_exp + '_mesh_hgr.nc'
144                 ENDIF ELSE BEGIN
145                    msg = 'eee : unset DRAKKAR_EXP keyword'
146                    ras = report(msg)
147                    msg = 'eee : orcares must be G42 or G70'
148                    ras = report(msg)
149                    RETURN
150                 ENDELSE
151              END
152      ELSE  : BEGIN
153                 msg = 'eee : invalid orcares parameter = ' + orcares
154                 ras = report(msg)
155                 msg = 'eee : orcares must be ORCA2 or ORCA025'
156                 ras = report(msg)
157                 RETURN
158              END
159  ENDCASE
160;
161; check method definition
162  CASE method OF
163     'bilinear': BEGIN
164                    msg = 'iii : valid method parameter = ' + method
165                    ras = report(msg)
166                 END
167     ELSE      : BEGIN
168                    msg = 'eee : invalid method parameter = ' + method
169                    ras = report(msg)
170                    msg = 'eee : method must be bilinear'
171                    ras = report(msg)
172                    RETURN
173                 END
174  ENDCASE
175;
176; check gridtype definition
177  gridtype = strupcase(gridtype)
178  CASE gridtype OF
179     'T' : BEGIN
180              msg = 'iii : valid gridtype parameter = ' + gridtype
181              ras = report(msg)
182           END
183     'U' : BEGIN
184              msg = 'iii : valid gridtype parameter = ' + gridtype
185              ras = report(msg)
186           END
187     'V' : BEGIN
188              msg = 'iii : valid gridtype parameter = ' + gridtype
189              ras = report(msg)
190           END
191     ELSE: BEGIN
192              msg = 'eee : invalid gridtype parameter = ' + gridtype
193              ras = report(msg)
194              msg = 'eee : gridtype must be T, U or V'
195              ras = report(msg)
196              RETURN
197           END
198  ENDCASE
199;
200; check for input files
201;
202; test if ${GEOMAG_ID} defined
203  geomag_id_env=GETENV('GEOMAG_ID')
204  CASE geomag_id_env OF
205     ''  :  BEGIN
206              msg = 'eee : ${GEOMAG_ID} is not defined'
207              ras = report(msg)
208              RETURN
209            END
210     ELSE: BEGIN
211             msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env
212             ras = report(msg)
213           END
214  ENDCASE
215;
216  iodirin = isadirectory(geomag_id_env)
217;
218; existence and protection of ${GEOMAG_ID}
219  IF (FILE_TEST(iodirin, /DIRECTORY,/EXECUTABLE , /READ) EQ 0) THEN BEGIN
220     msg = 'eee : the directory' + iodirin  + ' is not accessible.'
221     ras = report(msg)
222     RETURN
223  ENDIF
224;
225; build input filenames
226  filename_cond_sed = 'cond_sed' + '_' + orcares +'.nc'
227  fullfilename_cond_sed = iodirin + filename_cond_sed
228;
229; check if this file exists
230  fullfilename_cond_sed = isafile(iodirin + filename_cond_sed, NEW=0,/MUST_EXIST)
231  IF fullfilename_cond_sed[0] EQ '' THEN BEGIN
232     msg = 'eee : the file ' + fullfilename_cond_sed + ' was not found.'
233     ras = report(msg)
234     RETURN
235  ENDIF
236;
237; protection
238  IF (FILE_TEST(fullfilename_cond_sed[0], /READ) EQ 0) THEN BEGIN
239     msg = 'eee : the file ' + fullfilename_cond_sed[0] + ' is not readable.'
240     ras = report(msg)
241     RETURN
242  ENDIF
243;
244  filename_Br = 'Br' + '_' + orcares +'.nc'
245  fullfilename_Br = iodirin + filename_Br
246;
247; check if this file exists
248  fullfilename_Br = isafile(iodirin + filename_Br, NEW=0,/MUST_EXIST)
249  IF fullfilename_Br[0] EQ '' THEN BEGIN
250     msg = 'eee : the file ' + fullfilename_Br + ' was not found.'
251     ras = report(msg)
252     RETURN
253  ENDIF
254;
255; protection
256  IF (FILE_TEST(fullfilename_Br[0], /READ) EQ 0) THEN BEGIN
257     msg = 'eee : the file ' + fullfilename_Br[0] + ' is not readable.'
258     ras = report(msg)
259     RETURN
260  ENDIF
261;
262; mesh mask
263; check if this file exists
264  fullfilename_oce = isafile(iodirin + filename_oce, NEW=0, /MUST_EXIST, $
265                     RECURSIVE=0)
266  IF fullfilename_oce[0] EQ '' THEN BEGIN
267     msg = 'eee : the file ' + fullfilename_oce + ' was not found.'
268     ras = report(msg)
269     RETURN
270  ENDIF
271;
272; protection
273  IF (FILE_TEST(fullfilename_oce[0], /READ) EQ 0) THEN BEGIN
274     msg = 'eee : the file ' + fullfilename_oce[0] + ' is not readable.'
275     ras = report(msg)
276     RETURN
277  ENDIF
278;
279  filename_condmag = 'condmag.nc'
280  fullfilename_condmag = iodirin + filename_condmag
281;
282; check if this file exists
283  fullfilename_condmag = isafile(iodirin + filename_condmag, NEW=0,/MUST_EXIST)
284  IF fullfilename_condmag[0] EQ '' THEN BEGIN
285     msg = 'eee : the file ' + fullfilename_condmag + ' was not found.'
286     ras = report(msg)
287     RETURN
288  ENDIF
289;
290; protection
291  IF (FILE_TEST(fullfilename_condmag[0], /READ) EQ 0) THEN BEGIN
292     msg = 'eee : the file ' + fullfilename_condmag[0] + ' is not readable.'
293     ras = report(msg)
294     RETURN
295  ENDIF
296;
297; test if ${GEOMAG_OD} defined
298  geomag_od_env=GETENV('GEOMAG_OD')
299  CASE geomag_od_env OF
300     '' : BEGIN
301             msg = 'eee : ${GEOMAG_OD} is not defined'
302             ras = report(msg)
303             RETURN
304          END
305     ELSE: BEGIN
306             msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env
307             ras = report(msg)
308           END
309  ENDCASE
310;
311; check if output data will be possible
312  iodirout = isadirectory(geomag_od_env)
313;
314; existence and protection
315  IF (FILE_TEST(iodirout, /DIRECTORY,/WRITE) EQ 0) THEN BEGIN
316     msg = 'eee : the directory' + iodirout  + ' was not found.'
317     ras = report(msg)
318     RETURN
319  ENDIF
320;
321; build output filename
322  filename = 'condmag_from_' +  orcares +'.nc'
323  fullfilename = iodirout + filename
324;
325; in order to avoid unexpected overwritten
326  IF (FILE_TEST(fullfilename) EQ 1) THEN BEGIN
327     msg = 'eee : the file ' + fullfilename  + ' already exists.'
328     ras = report(msg)
329     RETURN
330  ENDIF
331;
332;----
333; Oceanic grid parameters ie input grid + mask
334;----
335;
336initocemesh, orcares, DRAKKAR_EXP = drakkar_exp
337;
338;----
339; read data to interpolate on regular grid
340;----
341;
342a=read_ncdf('cond_sed', 0, /timestep, $
343file = fullfilename_cond_sed[0], /nostruct)
344msg = 'iii : ' + fullfilename_cond_sed[0] + ' opened for read'
345ras = report(msg)
346; mask and limite cond_sed on values below 5000 ++
347cond_sedin = extrapsmooth(a, tmask[*, *, 0]*(a le 5000), /x_periodic)
348;
349b = read_ncdf('Br', 0, /timestep, $
350file = fullfilename_Br[0], /nostruct)
351msg = 'iii : ' + fullfilename_Br[0] + ' opened for read'
352ras = report(msg)
353; mask and limit Br on values below/above ++
354Brin = extrapsmooth(b, tmask[*, *, 0]*(b), /x_periodic)
355;
356glamin = glamt
357gphiin = gphit
358maskin = tmask[*, *, 0]
359;
360;---
361; conductivity and magnetic field grid parameters ie output grid
362;----
363
364initncdf, fullfilename_condmag[0], xaxisname = 'lo', yaxisname = 'la'
365msg = 'iii : ' + fullfilename_condmag[0] + ' opened for read'
366ras = report(msg)
367;
368glamout = glamt
369gphiout = gphit
370;
371
372dimidxout = jpi
373dimidyout = jpj
374;
375;---------------
376; Interpolate
377;---------------
378;++ cond_sedout = fromirr('bilinear', cond_sedin, glamin, gphiin, maskin,  glamout, gphiout, -1,  WEIG = weights, ADDR = addresses)
379cond_sedout=dindgen(dimidxout,dimidyout)
380;Brout=fromirr('bilinear', Brin, weights, addresses) ; ++ pb doc % Variable is undefined: LONOUT
381;++Brout = fromirr('bilinear', Brin, glamin, gphiin, maskin,  glamout, gphiout, -1,  WEIG = weights, ADDR = addresses)
382;
383Brout=dindgen(dimidxout,dimidyout)
384
385;
386;---------------
387; Produce outputs
388;---------------
389;
390  netcdf_id = NCDF_CREATE(fullfilename, /clobber)
391  NCDF_CONTROL, netcdf_id, /NOFILL
392;
393; dimension
394  dimidxout = NCDF_DIMDEF(netcdf_id, 'la' ,  jpi)
395  dimidyout = NCDF_DIMDEF(netcdf_id, 'lo' ,  jpj)
396;
397; global attributes
398  NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite
399;++  NCDF_ATTPUT, netcdf_id, 'file_name'  , filename, /GLOBAL
400;++  NCDF_ATTPUT, netcdf_id, 'Title'      , title, /GLOBAL
401;
402; declaration of variables
403  varid = lonarr(4)
404;
405  varid[0] = NCDF_VARDEF(netcdf_id, 'la'  , [dimidxout], /FLOAT)
406  varid[1] = NCDF_VARDEF(netcdf_id, 'lo'  , [dimidyout], /FLOAT)
407  varid[2] = NCDF_VARDEF(netcdf_id, 'cond_sed'  , [dimidxout, dimidyout], /FLOAT)
408  NCDF_ATTPUT, netcdf_id, varid[2], 'units'    , 'siemens'
409  varid[3] = NCDF_VARDEF(netcdf_id, 'Br'  , [dimidxout, dimidyout], /FLOAT)
410  NCDF_ATTPUT, netcdf_id, varid[3], 'units'    , 'siemens'
411;
412  NCDF_CONTROL, netcdf_id, /ENDEF
413;
414  NCDF_VARPUT, netcdf_id, varid[0], dindgen(dimidxout,dimidyout)
415  NCDF_VARPUT, netcdf_id, varid[1], dindgen(dimidyout,dimidyout)
416;
417; write the data
418;
419  NCDF_VARPUT, netcdf_id, varid[2], cond_sedout
420;++  NCDF_VARPUT, netcdf_id, 'Br', Brout
421;---------------------------
422; close the netcdf files
423  NCDF_CLOSE, netcdf_id
424;
425  msg = 'iii : ' + fullfilename + ' created'
426  ras = report(msg)
427;
428IF keyword_set(perf) EQ 1 THEN BEGIN
429  msg = 'iii : report profiler results'
430  ras = report(msg)
431; ++ tri par ordre alpha , par ordre de pourcentage croissant ou décroissant
432; +++ d'utilisation
433  profiler,/REPORT
434;  shut down all profiling (according to
435; http://www.dfanning.com/code_tips/whyslow.html)
436  profiler, /CLEAR, /SYSTEM
437  profiler, /CLEAR, /RESET
438ENDIF
439;
440END
Note: See TracBrowser for help on using the repository browser.