source: trunk/src/interp_erai_ws.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 Id
File size: 13.5 KB
Line 
1;+
2;
3;
4; ==================
5; interp_erai_ws.pro
6; ==================
7;
8; .. function:: interp_erai_ws(yyyymmddb,yyyymmdde, eraitype)
9;
10; DESCRIPTION
11; ===========
12;
13; Interpolation of u10 and v10 from ERA-I grid to OAFLUX grid
14;
15; If eraitype is set to 1,
16; :file:`${PROJECT_ID}/20c3m_erai_u10_TROP_1989_2009.nc`
17; containing
18; u10 from ERA-I
19; has been produced by
20; :ref:`compute_erai_daily_region_2d.sh`.
21;
22; If eraitype is set to 1,
23; :file:`${PROJECT_ID}/20c3m_erai_v10_TROP_1989_2009.nc`
24; containing
25; v10 from ERA-I
26; has been produced by
27; :ref:`compute_erai_daily_region_2d.sh`.
28;
29; If eraitype is set to 2,
30; :file:`${PROJECT_ID}/erai_u10_{yyyyb}_{yyyye}.nc`
31; containing
32; u10 from ERA-I
33; has been produced by
34; :ref:`concat_eraiy.sh`.
35;
36; If eraitype is set to 2,
37; :file:`${PROJECT_ID}/erai_v10_{yyyyb}_{yyyye}.nc`
38; containing
39; v10 from ERA-I
40; has been produced by
41; :ref:`concat_eraiy.sh`.
42;
43; :file:`${PROJECT_ID}/mask_oaflux_30N30S.nc`
44; containing OAFLUX grid
45; has been produced by
46; :func:`oaflux_mask_30n30s`.
47;
48; Interpolated u10 and v10 are written in
49; :file:`${PROJECT_OD}/erai_ws_19890101_20091231_oafluxgrid.nc`
50; if this file not already exists.
51;
52; This output file
53; :file:`${PROJECT_OD}/erai_ws_19890101_20091231_oafluxgrid.nc`
54; must be processed after by
55; :func:`ws_correction_ncdf`.
56;
57; .. only:: man
58;
59;    Figure is visible on PDF and HTML documents only.
60;
61; .. only:: html or latex
62;
63;     .. graphviz::
64;
65;        digraph interp_erai_ws {
66;
67;           file_u10_erai1 [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/20c3m_erai_u10_TROP_1989_2009.nc"];
68;           file_v10_erai1 [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/20c3m_erai_v10_TROP_1989_2009.nc"];
69;           file_u10_erai2 [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/erai_u10_{yyyyb}_{yyyye}.nc"];
70;           file_v10_erai2 [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/erai_v10_{yyyyb}_{yyyye}.nc"];
71;           mask [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/mask_oaflux_30N30S.nc"];
72;           file_out [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/erai_ws_19890101_20091231_oafluxgrid.nc"];
73;
74;           interp_erai_ws [shape=box,
75;           fontname=Courier,
76;           color=blue,
77;           URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/interp_erai_ws.pro",
78;           label="${PROJECT}/src/interp_erai_ws.pro"];
79;
80;           {file_u10_erai1 file_v10_erai1 file_u10_erai2 file_v10_erai2 mask} -> {interp_erai_ws} -> {file_out}
81;
82;        }
83;
84; SEE ALSO
85; ========
86;
87; :ref:`interpolate_data`
88;
89; :ref:`project_profile.sh`
90;
91; Used by :ref:`tropflux.sh`
92;
93; :ref:`compute_erai_daily_region_2d.sh`
94; :ref:`concat_eraiy.sh`
95;
96; :func:`report <saxo:report>`
97; :func:`isadirectory <saxo:isadirectory>`
98; :func:`isafile <saxo:isafile>`
99; :func:`initncdf <saxo:initncdf>`
100; :func:`read_ncdf <saxo:read_ncdf>`
101; :func:`domdef <saxo:domdef>`
102; :func:`call_interp2d <saxo:call_interp2d>`
103; :func:`ncdf_quickwrite <saxo:ncdf_quickwrite>`
104;
105; :func:`ws_correction_ncdf`
106;
107; EXAMPLES
108; ========
109;
110; With ERA-I type 1:
111;
112; .. code-block:: idl
113;
114;    .compile file_interp
115;    yyyymmddb = 20000101L
116;    yyyymmdde = 20001231L
117;    eraitype = 1
118;    result = interp_erai_ws(yyyymmddb, yyyymmdde, eraitype)
119;    print, result
120;
121; With ERA-I type 2:
122;
123; .. code-block:: idl
124;
125;    .compile file_interp
126;    yyyymmddb = 20000101L
127;    yyyymmdde = 20001231L
128;    eraitype = 2
129;    result = interp_erai_ws(yyyymmddb, yyyymmdde, eraitype)
130;
131; TODO
132; ====
133;
134; compare time between u10 and v10 after reading
135;
136; make it work : pb on loholt1 idl8
137;
138; reach end but bad values et nb timestep 4581:
139;
140; .. code-block:: bash
141;
142;    ncks -v u10  -d time,0,1 -d latitude,0,1 -d longitude,0,1 /homedata/pinsard/tropflux_d/erai_ws_19890101_20091231_oafluxgrid.nc
143;    latitude[0]=-29.5 degrees_north
144;    latitude[1]=-28.5 degrees_north
145;
146;    longitude[0]=30.5 degrees_east
147;    longitude[1]=31.5 degrees_east
148;
149;    time[0] latitude[0]=-29.5 longitude[0]=30.5 u10[0]=9.96921e+36 m/s
150;    time[0] latitude[0]=-29.5 longitude[1]=31.5 u10[1]=9.96921e+36 m/s
151;    time[0] latitude[1]=-28.5 longitude[0]=30.5 u10[350]=9.96921e+36 m/s
152;    time[0] latitude[1]=-28.5 longitude[1]=31.5 u10[351]=9.96921e+36 m/s
153;    time[1] latitude[0]=-29.5 longitude[0]=30.5 u10[21000]=9.96921e+36 m/s
154;    time[1] latitude[0]=-29.5 longitude[1]=31.5 u10[21001]=9.96921e+36 m/s
155;    time[1] latitude[1]=-28.5 longitude[0]=30.5 u10[21350]=9.96921e+36 m/s
156;    time[1] latitude[1]=-28.5 longitude[1]=31.5 u10[21351]=9.96921e+36 m/s
157;
158; make it work : pb on loholt1 idl6:
159;
160; .. parsed-literal::
161;
162;    Writing fields:
163;    u10[longitude,latitude,time]=u10out:u10_attr; v10[longitude,latitude,time]=v10out:v10_attr; longitude[]=xlon:lon_attr; latitude[]=ylat:lat_attr; tt[*time]=tt:time_attr  @ globattr
164;    % NCDF_VARPUT: Operation Failed, bad file (6) or variable [0] id ?
165;                   (NC_ERROR=-31)
166;    % Stop encountered: INTERP_ERAI_WS_1989_2009    1
167;
168; check OUTMASK_IND and SET_OUTMSKVAL added to call_interp2d call
169;
170; use real output of :ref:`compute_erai_daily_region_2d.sh`.
171;
172; coding rules
173;
174; check [yyyymmddb,yyyymmdde] validity
175;
176; check if [yyyymmddb,yyyymmdde] is included in ERA-I file
177;
178; what happen if yyyymmdde > 20091231 : need to read an other ERA-I file
179;
180; KNOWN ISSUES
181; ============
182;
183; test of existence of fullfilename_msk not very efficient because
184; MUST_EXIST keyword of :func:`isafile <saxo:isafile>` not yet implemented
185;
186; EVOLUTIONS
187; ==========
188;
189; $Id$
190;
191; $URL: svn+ssh://pinsard@forge.ipsl.jussieu.fr/ipsl/forge/projets/tropflux/svn/trunk/src/interp_erai_ws.pro $
192;
193; - fplod 20120704T113504Z cratos (Linux)
194;
195;   * add eraitype parameter
196;   * deal with time units (days vs hours) and time origin
197;
198; - fplod 20120320
199;
200;   * taking project_overwrite into account
201;   * try to add compile_opt seems to be incompatible with ncdf_quickwrite
202;   * pro -> function
203;   * hard coded st and en replaced by yyyymmddb and yyyymmdde parameters
204;   * start to homogenize time handling regarding :func:`interp_erai_t2m`
205;
206; - fplod 20110830T135212Z cratos (Linux)
207;
208;   * replace tt by time
209;
210; - pinsard 2011-08-23T09:09:16Z loholt1.ipsl.polytechnique.fr (Linux)
211;
212;   * retry on loholt1 with idl8 and idl6
213;
214; - fplod 20110808T094156Z cratos (Linux)
215;
216;   * add OUTMASK_IND and SET_OUTMSKVAL to call_interp2d call
217;
218; - pinsard 2011-07-05T07:33:36Z loholt1.ipsl.polytechnique.fr (Linux)
219;
220;   * usage of $PROJECT_ID and $PROJECT_OD
221;
222; - fplod 20101217T140745Z aedon.locean-ipsl.upmc.fr (Darwin)
223;
224;   * remove hard coded directory for mask_oaflux_30N30S.nc
225;
226; - fplod 20101215T112657Z aedon.locean-ipsl.upmc.fr (Darwin)
227;
228;   * add graph in header
229;
230; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin)
231;
232;   * minimal header
233;
234; - pbk 2008
235;
236;   * creation
237;
238;-
239function interp_erai_ws $
240         , yyyymmddb $
241         , yyyymmdde $
242         , eraitype
243;
244;++compile_opt idl2, strictarrsubs, logical_predicate
245;
246; Return to caller if errors
247ON_ERROR, 2
248;
249result = -1
250;
251usage = 'result = interp_erai_lwr(yyyymmddb, yyyymmdde, eraitype)'
252nparam = N_PARAMS()
253IF (nparam NE 3) THEN BEGIN
254    ras = report(['Incorrect number of arguments.' $
255    + '!C' $
256    + 'Usage : ' + usage])
257    return, result
258ENDIF
259;
260@cm_4cal
261@cm_4data
262@cm_4mesh
263@cm_4data
264@cm_project
265;
266; check for input directory
267;
268; test if ${PROJECT_ID} defined
269CASE project_id_env OF
270    '' : BEGIN
271        msg = 'eee : ${PROJECT_ID} is not defined'
272        ras = report(msg)
273        return, result
274    END
275    ELSE : BEGIN
276        msg = 'iii : ${PROJECT_ID} is ' + project_id_env
277        ras = report(msg)
278    END
279ENDCASE
280;
281iodirin = isadirectory(project_id_env)
282;
283; existence and protection of ${PROJECT_ID}
284IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
285    msg = 'eee : the directory' + iodirin  + ' is not accessible.'
286    ras = report(msg)
287    return, result
288ENDIF
289;
290; build mask filename
291filename_msk='mask_oaflux_30N30S.nc'
292;
293; check if this file exists
294msg='iii : looking for ' + filename_msk
295ras = report(msg)
296fullfilename_msk = isafile(iodirin + filename_msk, NEW=0, /MUST_EXIST)
297IF fullfilename_msk[0] EQ '' THEN BEGIN
298    msg = 'eee : the file ' + fullfilename_msk + ' was not found.'
299    ras = report(msg)
300    return, result
301ENDIF
302;
303; build u10 data filename
304IF eraitype EQ 1 THEN BEGIN
305    erai_var = 'u10'
306    filename='20c3m_erai_u10_TROP_1989_2009.nc'
307    tutoday = 24.d
308    tmporig = date2jul(19570101L,month=immo,day=iddo,year=iyyyyo)
309endif
310IF (eraitype eq 2) THEN BEGIN
311    erai_var = 'u10'
312    tmpb = date2jul(yyyymmddb,month=immb,day=iddb,year=iyyyyb)
313    tmpe = date2jul(yyyymmdde,month=imme,day=idde,year=iyyyye)
314    filename = 'erai_' + erai_var + '_' + string(iyyyyb,format='(I4.4)') + '_' + string(iyyyye,format='(I4.4)') + '.nc'
315    tutoday = 1.d
316    tmporig = date2jul(19500101L,month=immo,day=iddo,year=iyyyyo)
317ENDIF
318;
319; check if this file exists
320msg='iii : looking for ' + filename
321ras = report(msg)
322fullfilename = isafile(iodirin + filename, NEW=0, /MUST_EXIST)
323IF fullfilename[0] EQ '' THEN BEGIN
324    msg = 'eee : the file ' + fullfilename + ' was not found.'
325    ras = report(msg)
326    return, result
327ENDIF
328;
329; test if ${PROJECT_OD} defined
330CASE project_od_env OF
331    '' : BEGIN
332        msg = 'eee : ${PROJECT_OD} is not defined'
333        ras = report(msg)
334        return, result
335    END
336    ELSE : BEGIN
337        msg = 'iii : ${PROJECT_OD} is ' + project_od_env
338        ras = report(msg)
339    END
340ENDCASE
341;
342; check if output data will be possible
343iodirout = isadirectory(project_od_env)
344;
345; existence and protection
346IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN
347    msg = 'eee : the directory' + iodirout  + ' was not found.'
348    ras = report(msg)
349    return, result
350ENDIF
351;
352; build output filename
353filename_out = 'erai_ws_' + string(yyyymmddb,format='(I8.8)') + '_' + string(yyyymmdde,format='(I8.8)') + '_oafluxgrid.nc'
354fullfilename_out = iodirout + filename_out
355; in order to avoid unexpected overwritten
356IF ((FILE_TEST(fullfilename_out) EQ 1)  AND (project_overwrite EQ 0)) THEN BEGIN
357    msg = 'eee : the file ' + fullfilename_out  + ' already exists.'
358    ras = report(msg)
359    return, result
360ENDIF
361;
362initncdf, fullfilename
363u10in=read_ncdf(erai_var, yyyymmddb - .5d, yyyymmdde, file=fullfilename,/nostr)
364;
365timein=tutoday*(time-julday(immo,iddo,iyyyyo,0,0,0))
366jpt=n_elements(timein)
367da=jul2date(time[0])
368cda0=string(da,format='(i8.8)')
369da=jul2date(time[jpt-1])
370cda1=string(da,format='(i8.8)')
371print, 'u10 first date ', cda0
372print, 'u10 last date ' , cda1
373;
374; build v10 data filename
375IF eraitype EQ 1 THEN BEGIN
376    erai_var = 'v10'
377    filename='20c3m_erai_v10_TROP_1989_2009.nc'
378endif
379IF (eraitype eq 2) THEN BEGIN
380    erai_var = 'v10'
381    tmpb = date2jul(yyyymmddb,month=immb,day=iddb,year=iyyyyb)
382    tmpe = date2jul(yyyymmdde,month=imme,day=idde,year=iyyyye)
383    filename = 'erai_' + erai_var + '_' + string(iyyyyb,format='(I4.4)') + '_' + string(iyyyye,format='(I4.4)') + '.nc'
384ENDIF
385;
386; check if this file exists
387msg='iii : looking for ' + filename
388ras = report(msg)
389fullfilename = isafile(iodirin + filename, NEW=0, /MUST_EXIST)
390IF fullfilename[0] EQ '' THEN BEGIN
391    msg = 'eee : the file ' + fullfilename + ' was not found.'
392    ras = report(msg)
393    return, result
394ENDIF
395initncdf, fullfilename
396v10in=read_ncdf(erai_var, yyyymmddb - .5d, yyyymmdde, file=fullfilename,/nostr)
397;
398initncdf, fullfilename
399domdef
400latin=reform(gphit[0,*])
401lonin=reform(glamt[*,0])
402print, 'lat grid ',min(latin),max(latin),latin[1]-latin[0]
403print, 'lon grid ',min(lonin),max(lonin),lonin[1]-lonin[0]
404;
405timein=tutoday*(time-julday(immo,iddo,iyyyyo,0,0,0))
406jptin=n_elements(timein)
407da=jul2date(time[0])
408cda0=string(da,format='(i8.8)')
409da=jul2date(time[jpt-1])
410cda1=string(da,format='(i8.8)')
411print, 'v10 first date ', cda0
412print, 'v10 last date ' , cda1
413help, v10in
414;
415tab=u10in[*,*,0]
416mskin=glamt*0.+1.
417;
418initncdf, fullfilename_msk
419domdef
420latout=reform(gphit[0,*])
421lonout=reform(glamt[*,0])
422print, 'lat grid ',min(latout),max(latout),latout[1]-latout[0]
423print, 'lon grid ',min(lonout),max(lonout),lonout[1]-lonout[0]
424mskout=read_ncdf("msk", file=fullfilename_msk,/nostr)
425;
426help, u10in, v10in,lonin,latin,mskin,lonout,latout,mskout
427;
428u10out=fltarr(jpi,jpj,jptin)
429v10out=fltarr(jpi,jpj,jptin)
430;
431for jt=0,jptin-1 do begin
432    ; ++print, 'Interpolation jt=',jt,' / ',jptin-1
433    tab=reform(u10in[*,*,jt])
434    u10out[*,*,jt]=call_interp2d(tab,lonin,latin,mskin $
435    , lonout,latout,method='bilinear' $
436    , OUTMASK_IND=mskout, SET_OUTMSKVAL=mskout)
437    u10out[*,*,jt]=u10out[*,*,jt]*mskout+(1.-mskout)*1.e20
438    ;
439    tab=reform(v10in[*,*,jt])
440    v10out[*,*,jt]=call_interp2d(tab,lonin,latin,mskin $
441    , lonout,latout,method='bilinear' $
442    , OUTMASK_IND=mskout, SET_OUTMSKVAL=mskout)
443    v10out[*,*,jt]=v10out[*,*,jt]*mskout+(1.-mskout)*1.e20
444    ;
445endfor
446;
447xlon=reform(glamt[*,0] )
448ylat=reform(gphit[0,*])
449;
450initncdf, fullfilename_msk
451valmask=1.e20
452;
453ylat=latout
454xlon=lonout
455ncfile='!' + fullfilename_out
456lon_attr={units:'degrees_east',long_name:'Longitude'}
457lat_attr={units:'degrees_north',long_name:'Latitude'}
458time_attr={units:'days since 1950-01-01 00:00:00',long_name:'Time axis',time_origin:'1950-JAN-01 00:00:00'}
459globattr={source:'Tropical ocean winds obtained from ERA Interim',timerange:cda0+' - '+cda1}
460u10_attr={units:'m/s',missing_value:valmask,long_name:'10 metre u wind component',short_name:'u10',axis:'TYX'}
461v10_attr={units:'m/s',missing_value:valmask,long_name:'10 metre v wind component',short_name:'v10',axis:'TYX'}
462;
463help, u10out
464help, v10out
465help, timein
466;
467ncfields = 'u10[longitude,latitude,*time]=u10out:u10_attr; ' $
468+ 'v10[longitude,latitude,*time]=v10out:v10_attr; ' $
469+ 'longitude[]=xlon:lon_attr; ' $
470+ 'latitude[]=ylat:lat_attr; ' $
471+ 'time[]=timein:time_attr ' $
472+ ' @ globattr'
473;
474@ncdf_quickwrite
475;
476result = 0
477return, result
478;
479end
Note: See TracBrowser for help on using the repository browser.