source: trunk/src/ws_correction_ncdf.pro

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

fix thanks to coding rules; typo

  • Property svn:keywords set to URL
File size: 7.1 KB
Line 
1;+
2;
3; ======================
4; ws_correction_ncdf.pro
5; ======================
6;
7; .. function:: ws_correction_ncdf(yyyymmddb,yyyymmdde)
8;
9; DESCRIPTION
10; ===========
11;
12; Correction of ws on OAFLUX grid
13;
14; :file:`${PROJECT_OD}/erai_ws_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc`
15; containing
16; ++
17; have been produced by
18; :func:`interp_erai_ws`.
19;
20; Corrected ws on OAFLUX grid
21; is written in
22; :file:`${PROJECT_OD}/TropFlux_ws_{yyyymmdd}_{yyyymmdd}.nc`
23; if this file not already exists.
24;
25; This output file
26; :file:`${PROJECT_OD}/TropFlux_ws_{yyyymmdd}_{yyyymmdd}.nc`
27; will be used by
28; :func:`tropflux`.
29;
30; .. only:: man
31;
32;    Figure is visible on PDF and HTML documents only.
33;
34; .. only:: html or latex
35;
36;     .. graphviz::
37;
38;        digraph ws_correction_ncdf {
39;
40;           file_in [shape=ellipse,
41;           fontname=Courier,
42;           label="${PROJECT_OD}/erai_ws_{yyyymmdd}_{yyyymmdd}_oafluxgrid.nc"];
43;           file_out [shape=ellipse,
44;           fontname=Courier,
45;           label="${PROJECT_OD}/TropFlux_ws_{yyyymmdd}_{yyyymmdd}.nc"];
46;
47;           ws_correction_ncdf [shape=box,
48;           fontname=Courier,
49;           color=blue,
50;           URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/ws_correction_ncdf.pro",
51;           label="${PROJECT}/src/ws_correction_ncdf.pro"];
52;
53;           {file_in} -> {ws_correction_ncdf} -> {file_out}
54;
55;        }
56;
57; SEE ALSO
58; ========
59;
60; :ref:`mooring_corrections`
61;
62; Used by :ref:`tropflux.sh`
63;
64; :ref:`project_profile.sh`
65;
66; :func:`interp_erai_ws`
67;
68; :func:`initncdf <saxo:initncdf>`
69; :func:`read_ncdf <saxo:read_ncdf>`
70; :func:`grossemoyenne <saxo:grossemoyenne>`
71; :func:`julday <saxo:julday>`
72; :func:`caldat <saxo:caldat>`
73; :func:`ncdf_quickwrite <saxo:ncdf_quickwrite>`
74;
75; EXAMPLES
76; ========
77;
78; .. code-block:: idl
79;
80;    yyyymmddb = 20000101L
81;    yyyymmdde = 20001231L
82;    result = ws_correction_ncdf(yyyymmddb, yyyymmdde)
83;    print, result
84;
85; TODO
86; ====
87;
88; coding rules
89;
90; hard coded correction
91;
92; KNOWN ISSUES
93; ============
94;
95; test of existence of fullfilename_in not very efficient because
96; MUST_EXIST keyword of :func:`isafile <saxo:isafile>` not yet implemented
97;
98; EVOLUTIONS
99; ==========
100;
101; $Id: ws_correction_ncdf.pro 88 2011-08-19 15:40:14Z pinsard $
102;
103; $URL$
104;
105; - fplod 20120322
106;
107;   * taking project_overwrite into account
108;   * try to add compile_opt seems to be incompatible with ncdf_quickwrite
109;   * pro -> function
110;   * hard coded da1 and da2 replaced by yyyymmddb and yyyymmdde parameters
111;   * get rid of timegen
112;
113; - fplod 20110808T131954Z aedon.locean-ipsl.upmc.fr (Darwin)
114;
115;   * usage of ${PROJECT_OD}
116;   * complete description
117;   * remove v20 in output filename
118;
119; - fplod 20101215T115916Z aedon.locean-ipsl.upmc.fr (Darwin)
120;
121;   * add graph in header
122;
123; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin)
124;
125;   * minimal header
126;
127; - pbk 2008
128;
129;   * creation
130;
131;-
132function ws_correction_ncdf $
133, yyyymmddb $
134, yyyymmdde
135;
136;++compile_opt idl2, strictarrsubs, logical_predicate
137;
138@cm_4data
139@cm_4mesh
140@cm_4data
141@cm_project
142;
143; Return to caller if errors
144ON_ERROR, 2
145;
146result = -1
147;
148usage = 'result = ws_correction_ncdf(yyyymmddb, yyyymmdde)'
149nparam = N_PARAMS()
150IF (nparam NE 2) THEN BEGIN
151    ras = report(['Incorrect number of arguments.' $
152          + '!C' $
153          + 'Usage : ' + usage])
154    return, result
155ENDIF
156;
157@cm_4cal
158; test if ${PROJECT_OD} defined
159CASE project_od_env OF
160    '' : BEGIN
161        msg = 'eee : ${PROJECT_OD} is not defined'
162        ras = report(msg)
163        return, result
164    END
165    ELSE : BEGIN
166        msg = 'iii : ${PROJECT_OD} is ' + project_od_env
167        ras = report(msg)
168    END
169ENDCASE
170;
171; check if output data will be possible
172iodirout = isadirectory(project_od_env)
173;
174; existence and protection for reading
175IF (FILE_TEST(iodirout, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
176    msg = 'eee : the directory' + iodirout  + ' is not accessible.'
177    ras = report(msg)
178    return, result
179ENDIF
180;
181; existence and protection for writing
182IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN
183    msg = 'eee : the directory' + iodirout  + ' was not found.'
184    ras = report(msg)
185    return, result
186ENDIF
187;
188; build data filename
189filename='erai_ws_' $
190    + string(yyyymmddb,format='(i8.8)') $
191    + '_' $
192    + string(yyyymmdde,format='(i8.8)') $
193    + '_oafluxgrid.nc'
194;
195; check if this file exists
196msg='iii : looking for ' + filename
197ras = report(msg)
198fullfilename = isafile(iodirout + filename, NEW=0, /MUST_EXIST)
199IF fullfilename[0] EQ '' THEN BEGIN
200    msg = 'eee : the file ' + fullfilename + ' was not found.'
201    ras = report(msg)
202    return, result
203ENDIF
204;
205; build output filename
206filename_out = 'TropFlux_ws_' $
207    + string(yyyymmddb,format='(i8.8)') $
208    + '_' $
209    + string(yyyymmdde,format='(i8.8)') $
210    + '.nc'
211fullfilename_out = iodirout + filename_out
212; in order to avoid unexpected overwritten
213IF ((FILE_TEST(fullfilename_out) EQ 1)  AND (project_overwrite EQ 0)) THEN BEGIN
214    msg = 'eee : the file ' + fullfilename_out  + ' already exists.'
215    ras = report(msg)
216    return, result
217ENDIF
218;
219initncdf, fullfilename
220u=read_ncdf('u10',yyyymmddb-.5d,yyyymmdde,file=fullfilename,/nostr)
221timein=24.d*(time-julday(1,1,1957,0,0,0))
222jpt=n_elements(timein)
223da=jul2date(time[0])
224cda0=string(da,format='(i8.8)')
225da=jul2date(time[jpt-1])
226cda1=string(da,format='(i8.8)')
227print, 'u10 in ws_correction_ncdf first date ', cda0
228print, 'u10 in ws_correction_ncdf last date ' , cda1
229;
230v=read_ncdf('v10',yyyymmddb-.5d,yyyymmdde,file=fullfilename,/nostr)
231timein=24.d*(time-julday(1,1,1957,0,0,0))
232jpt=n_elements(timein)
233da=jul2date(time[0])
234cda0=string(da,format='(i8.8)')
235da=jul2date(time[jpt-1])
236cda1=string(da,format='(i8.8)')
237print, 'v10 in ws_correction_ncdf first date ', cda0
238print, 'v10 in ws_correction_ncdf last date ' , cda1
239;
240w=sqrt(u*u+v*v)
241;
242help, w
243;
244w_mean=grossemoyenne(w,'t',/nan)
245help, w_mean
246;
247jpt=n_elements(time)
248caldat, time,mon,day,yea
249w_m=w*0.
250;
251for jt=0,jpt-1 do begin
252    jtt=(time[jt]-julday(1,1,yea[jt])) < 364
253    q=reform(w_mean[*,*])
254    w_m[*,*,jt]=q
255endfor
256help, w_m
257;
258w_ano=w-w_m
259;
260; correction for mean based on scatter
261;w_m=w_m+0.282667    ; (2000-2008)
262w_m=w_m+0.276739     ; (2000-2009)
263;
264help, w_ano
265;
266; applying the correction for variability based on the scatter
267;w_ano=w_ano*(1/0.897667)   ; (2000-2008)
268w_ano=w_ano*(1/0.903587)    ; (2000-2009)
269;
270w_new=w_m+w_ano
271help, w_new
272;
273;writing field
274lat=reform(gphit[0,0:jpj-1])
275lon=reform(glamt[0:jpi-1,0])
276;
277ncfile='!' + fullfilename_out
278lon_attr={units:'degrees_east',long_name:'Longitude'}
279lat_attr={units:'degrees_north',long_name:'Latitude'}
280time_attr={units:'days since 1950-01-01 00:00:00',long_name:'Time axis',time_origin:'1950-JAN-01 00:00:00'}
281globattr={source:'Basic data obtained from ERAI. Bias and variability correction are applied',timerange:cda0+' - '+cda1}
282w_attr={units:'m/s',missing_value:1.e20,long_name:'mean wind speed',short_name:'w',axis:'TYX'}
283;
284ncfields = 'ws[longitude,latitude,*time]=w_new:w_attr; ' $
285                      + 'longitude[]=lon:lon_attr; ' $
286                      + 'latitude[]=lat:lat_attr; ' $
287                      + 'time[]=timein:time_attr ' $
288                      + ' @ globattr'
289;
290@ncdf_quickwrite
291;
292result = 0
293return, result
294;
295end
Note: See TracBrowser for help on using the repository browser.