source: trunk/step2_diff.pro @ 39

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

step2_diff ok
  • Property svn:keyword set to ID
File size: 7.1 KB
Line 
1;+
2; @file_comments
3; compute delta between *_5d_yyyy_grid_T_orcares.pro Netcdf files with same dimension
4;
5; @param file1 {in}{required} {type=string}
6; name of the first file to be read
7;
8; @param file2 {in}{required} {type=string}
9; name of the second file to be read
10;
11; @param file3 {in}{required} {type=string}
12; file where differences between variables are written
13;
14; @examples
15; test (for header and delta)
16; to compute difference between Sigma in Sigma_5d_1993_grid_T_ORCA2.nc and
17; itself and write it in ginette.nc :
18; $ cd ${GEOMAG}
19; $ idl
20; IDL> file1=getenv('GEOMAG_OD') + '/Sigma_5d_1993_grid_T_ORCA2.nc'
21; IDL> file2=getenv('GEOMAG_OD') + '/Sigma_5d_1993_grid_T_ORCA2.nc'
22; IDL> step2_diff, file1, file2, 'ginette.nc'
23; values of delta must be 0 everywhere
24; $ rm ginette.nc
25;
26; real life
27; to compute difference between Sigma in
28; /usr/work/sur/fvi/OPA/ORCA2/Sigma_5d_1994_grid_T.nc and
29; ${GEOMGA_OD}/Sigma_5d_1994_grid_T_ORCA2.nc and write it in
30; Sigma_5d_1994_grid_T_ORCA2_diff.nc :
31; $ cd ${GEOMAG}
32; $ idl
33; IDL> year=1994
34; IDL> file1= '/usr/work/sur/fvi/OPA/ORCA2/Sigma_5d_1994_grid_T.nc'
35; IDL> file2=getenv('GEOMAG_OD') + '/Sigma_5d_1994_grid_T_ORCA2.nc'
36; IDL> step2_diff, file1, file2, 'Sigma_5d_1994_grid_T_ORCA2_diff.nc'
37;
38; to see the difference, for example
39; IDL> xxx,'Sigma_5d_1994_grid_T_ORCA2_diff.nc'
40; select xy on plt wigdet
41;
42; see also step2_diff.sh for quick start
43;
44; @history
45; fplod 2007-08-22T09:39:45Z aedon.locean-ipsl.upmc.fr (Darwin)
46; correction
47; fplod 2007-07-30T10:58:12Z aedon.locean-ipsl.upmc.fr (Darwin)
48; like <progeomag>step1_diff</progeomag>
49;
50; @todo
51; check difference between headers
52;
53; @version
54; $Id$
55;
56PRO step2_diff, file1, file2, file3
57;
58ncverbose=1
59;
60; check if file1 exists
61fullfile1 = isafile(file1, NEW=0,/MUST_EXIST)
62IF fullfile1[0] EQ '' THEN BEGIN
63   msg = 'eee : the file ' + fullfile1 + ' was not found.'
64   ras = report(msg)
65   RETURN
66ENDIF
67;
68; protection
69IF (FILE_TEST(fullfile1[0], /READ) EQ 0) THEN BEGIN
70   msg = 'eee : the file ' + fullfile1[0] + ' is not readable.'
71   ras = report(msg)
72   RETURN
73ENDIF
74;
75; check if file2 exists
76fullfile2 = isafile(file2, NEW=0,/MUST_EXIST)
77IF fullfile2[0] EQ '' THEN BEGIN
78   msg = 'eee : the file ' + fullfile2 + ' was not found.'
79   ras = report(msg)
80   RETURN
81ENDIF
82;
83; protection
84IF (FILE_TEST(fullfile2[0], /READ) EQ 0) THEN BEGIN
85   msg = 'eee : the file ' + fullfile2[0] + ' is not readable.'
86   ras = report(msg)
87   RETURN
88ENDIF
89;
90; read the first file
91ncdf_read,  fullfile1[0], file1info, file1dinfo, file1vinfo, $
92            file1gatts, file1vatts, file1data
93;
94; read the second file
95ncdf_read,  fullfile2[0], file2info, file2dinfo, file2vinfo, $
96            file2gatts, file2vatts, file2data
97;
98; compute delta
99delta = *file1data[4] - *file2data[4]
100minarr = min(delta, max = maxarr)
101;
102; check for number of non zero in delta
103checknonzero=where(delta NE 0.,count)
104IF count EQ 0 THEN BEGIN
105   msg = 'iii : delta is zero everywhere'
106   ras = report(msg)
107ENDIF ELSE BEGIN
108   msg = 'iii : delta is not zero ' + STRING(count)  + ' times'
109   ras = report(msg)
110ENDELSE
111;
112; write output
113; in order to avoid unexpected overwritten
114IF (FILE_TEST(file3) EQ 1) THEN BEGIN
115   msg = 'eee : the file ' + file3  + ' already exists.'
116   ras = report(msg)
117   RETURN
118ENDIF
119;
120;---------------------------
121; Creation of the NetCdf output file
122;---------------------------
123;
124  netcdf_id = NCDF_CREATE(file3, /clobber)
125  NCDF_CONTROL, netcdf_id, /NOFILL
126;
127; dimension
128  dimidx = NCDF_DIMDEF(netcdf_id, file1dinfo[0].name ,  file1dinfo[0].size)
129  dimidy = NCDF_DIMDEF(netcdf_id, file1dinfo[1].name ,  file1dinfo[1].size)
130  dimidz = NCDF_DIMDEF(netcdf_id, file1dinfo[2].name ,  file1dinfo[2].size)
131  dimidt = NCDF_DIMDEF(netcdf_id, file1dinfo[3].name, /UNLIMITED)
132  dimidz2 = NCDF_DIMDEF(netcdf_id, file1dinfo[4].name ,  file1dinfo[4].size)
133;
134; global attributes
135  NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite
136  NCDF_ATTPUT, netcdf_id, 'file_name'  , file3, /GLOBAL
137  NCDF_ATTPUT, netcdf_id, 'Title'      , $
138               'delta between '+ file1 + ' and ' + file2, /GLOBAL
139;
140; declaration of variables
141  varid = lonarr(5)
142;
143  varid[0] = NCDF_VARDEF(netcdf_id, file1vinfo[0].name  , [dimidx, dimidy], /FLOAT)
144  NCDF_ATTPUT, netcdf_id, varid[0], $
145               (*file1vatts[0])[0].name , $
146               STRING(*(*file1vatts[0])[0].values)
147  NCDF_ATTPUT, netcdf_id, varid[0], $
148               (*file1vatts[0])[1].name, $
149               STRING(*(*file1vatts[0])[1].values)
150  NCDF_ATTPUT, netcdf_id, varid[0], $
151               (*file1vatts[0])[2].name, $
152               STRING(*(*file1vatts[0])[2].values)
153;
154  varid[1] = NCDF_VARDEF(netcdf_id, file1vinfo[1].name  , [dimidx, dimidy], /FLOAT)
155  NCDF_ATTPUT, netcdf_id, varid[1], $
156               (*file1vatts[1])[0].name, $
157               STRING(*(*file1vatts[1])[0].values)
158  NCDF_ATTPUT, netcdf_id, varid[1], $
159               (*file1vatts[1])[1].name, $
160               STRING(*(*file1vatts[1])[1].values)
161  NCDF_ATTPUT, netcdf_id, varid[1], $
162               (*file1vatts[1])[2].name, $
163               STRING(*(*file1vatts[1])[2].values)
164;
165  varid[2] = NCDF_VARDEF(netcdf_id, file1vinfo[2].name, [dimidz2], /FLOAT);
166  NCDF_ATTPUT, netcdf_id, varid[2], $
167               (*file1vatts[2])[0].name, $
168               STRING(*(*file1vatts[2])[0].values)
169  NCDF_ATTPUT, netcdf_id, varid[2], $
170               (*file1vatts[2])[1].name, $
171               STRING(*(*file1vatts[2])[1].values)
172  NCDF_ATTPUT, netcdf_id, varid[2], $
173               (*file1vatts[2])[2].name, $
174               STRING(*(*file1vatts[2])[2].values)
175;
176  varid[3] = NCDF_VARDEF(netcdf_id, file1vinfo[3].name, [dimidt], /FLOAT);
177  NCDF_ATTPUT, netcdf_id, varid[3], $
178               (*file1vatts[3])[0].name, $
179               STRING(*(*file1vatts[3])[0].values)
180  NCDF_ATTPUT, netcdf_id, varid[3], $
181               (*file1vatts[3])[1].name, $
182               STRING(*(*file1vatts[3])[1].values)
183  NCDF_ATTPUT, netcdf_id, varid[3], $
184               (*file1vatts[3])[2].name, $
185               STRING(*(*file1vatts[3])[2].values)
186  NCDF_ATTPUT, netcdf_id, varid[3], $
187               (*file1vatts[3])[3].name, $
188               STRING(*(*file1vatts[3])[3].values)
189  NCDF_ATTPUT, netcdf_id, varid[3], $
190               (*file1vatts[3])[4].name, $
191               STRING(*(*file1vatts[3])[4].values)
192help,delta
193
194  varid[4] =  NCDF_VARDEF(netcdf_id, 'delta',   [dimidx, dimidy,dimidz, dimidt], /DOUBLE)
195  NCDF_ATTPUT, netcdf_id, varid[4], 'long_name', 'delta'
196  NCDF_ATTPUT, netcdf_id, varid[4], 'valid_min', minarr ,/DOUBLE
197  NCDF_ATTPUT, netcdf_id, varid[4], 'valid_max', maxarr ,/DOUBLE
198;
199;---------------------------
200; end of header definition, writing of the NetCdf files
201;---------------------------
202;
203  NCDF_CONTROL, netcdf_id, /ENDEF
204;
205  NCDF_VARPUT, netcdf_id, file1vinfo[0].name, (*file1data[0])
206  NCDF_VARPUT, netcdf_id, file1vinfo[1].name, (*file1data[1])
207  NCDF_VARPUT, netcdf_id, file1vinfo[2].name, (*file1data[2])
208  NCDF_VARPUT, netcdf_id, file1vinfo[3].name, (*file1data[3])
209;
210; write the data
211  NCDF_VARPUT, netcdf_id, 'delta', delta
212;
213;---------------------------
214; close the netcdf files
215  NCDF_CLOSE, netcdf_id
216;
217  msg = 'iii : ' + file3 + ' created'
218  ras = report(msg)
219;
220END
Note: See TracBrowser for help on using the repository browser.