source: trunk/condmag_from_orca.pro

Last change on this file was 48, checked in by pinsard, 10 years ago

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