source: trunk/src/TropFlux_19890101_20091231.pro @ 79

Last change on this file since 79 was 77, checked in by pinsard, 13 years ago

complete documentation
  • Property svn:executable set to *
File size: 14.6 KB
Line 
1;+
2;
3; .. _TropFlux_19890101_20091231.pro:
4;
5; ==============================
6; TropFlux_19890101_20091231.pro
7; ==============================
8;
9; This program computes net heat flux components on the 1° oaflux grid.
10;
11; all input variables are corrected for mean bias and variability.
12;
13; gustiness correction is applied for wind speed based on cronin's climatological
14; gustiness values.
15;
16; :file:`${PROJECT_ID}/mask_oaflux_30N30S.nc` containing OAFLUX grid have been produced by :ref:`oaflux_mask_30N30S.pro`.
17;
18; :file:`${PROJECT_OD}/TropFlux_sst_19890101_20091231.nc` containing sst corrected on OAFLUX grid have been produced by :ref:`sst_correction_ncdf.pro`.
19;
20; :file:`${PROJECT_OD}/TropFlux_ws_19890101_20091231.nc` containing ws corrected on OAFLUX grid have been produced by :ref:`ws_correction_ncdf.pro`.
21;
22; :file:`${PROJECT_OD}/TropFlux_gustiness_19890101_20091231.nc` containing ++
23; have been produced by :ref:`cronin_gustiness_ncdf.pro`.
24;
25; :file:`${PROJECT_OD}/TropFlux_swr_19890101_20091231_BLND.nc` containing ws corrected on OAFLUX grid
26; have been produced by :ref:`TropFlux_swr_BLND_19890101_20091231.pro`.
27;
28; :file:`${PROJECT_OD}/TropFlux_lwr_19890101_20091231.nc` containing lwr corrected on OAFLUX grid have been produced by :ref:`lwr_correction_ncdf.pro`.
29;
30; :file:`${PROJECT_OD}/TropFlux_t2m_19890101_20091231.nc` containing t2m corrected on OAFLUX grid have been produced by :ref:`t2m_correction_ncdf.pro`.
31;
32; :file:`${PROJECT_OD}/TropFlux_q2m_19890101_20091231.nc` containing q2m corrected on OAFLUX grid have been produced by :ref:`d2m_to_q2m_erai.pro`.
33;
34; net heat flux components are written
35; in :file:`${PROJECT_OD}/TropFlux_19890101_20091231.nc`
36; if this file not already exists.
37;
38;     .. graphviz::
39;
40;        digraph tropflux_19890101_20091231 {
41;           graph [
42;           rankdir="LR",
43;           ]
44;
45;           mask [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/mask_oaflux_30N30S.nc"];
46;           file_sst [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_sst_19890101_20091231.nc"];
47;           file_ws [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_ws_19890101_20091231.nc"];
48;           file_wg [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_gustiness_19890101_20091231.nc"];
49;           file_swr [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_swr_19890101_20091231_BLND.nc"];
50;           file_lwr [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_lwr_19890101_20091231.nc"];
51;           file_t2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_t2m_19890101_20091231.nc"];
52;           file_q2m [shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_q2m_19890101_20091231.nc"];
53;
54;           file_out[shape=ellipse,fontname=Courier,label="${PROJECT_OD}/TropFlux_19890101_20091231.nc"];
55;
56;           tropflux_19890101_20091231 [shape=box,
57;           fontname=Courier,
58;           color=blue,
59;           URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/TropFlux_19890101_20091231.pro",
60;           label="${PROJECT}/src/TropFlux_19890101_20091231.pro"];
61;
62;           {mask file_sst file_ws file_wg file_swr file_lwr file_t2m file_q2m} -> {tropflux_19890101_20091231} -> {file_out}
63;
64;          }
65;
66; SEE ALSO
67; ========
68;
69; :ref:`project_profile.sh`
70;
71; :func:`report <saxo:report>`
72; :func:`isadirectory <saxo:isadirectory>`
73; :func:`isafile <saxo:isafile>`
74; :func:`initncdf <saxo:initncdf>`
75; :func:`ncdf_lec <saxo:ncdf_lec>`
76; :func:`read_ncdf <saxo:read_ncdf>`
77; :func:`caldat <saxo:caldat>`
78; :func:`julday <saxo:julday>`
79; :func:`ncdf_quickwrite <saxo:ncdf_quickwrite>`
80;
81; :func:`cor30a`
82;
83; EXAMPLES
84; ========
85;
86; ::
87;
88;  IDL> .compile TropFlux_19890101_20091231
89;  IDL> tropflux_19890101_20091231
90;
91; TODO
92; ====
93;
94; I (fp) do not know how
95; ${PROJECT_OD}/TropFlux_swr_19890101_20091231_BLND.nc is produced.
96;
97; avoid mix lower/uppercase in pro name to avoid compile
98;
99; coding rules
100;
101; KNOWN ISSUES
102; ============
103;
104; test of existence of fullfilename_msk not very efficient because
105; MUST_EXIST keyword of :func:`isafile <saxo:isafile>` not yet implemented
106;
107; EVOLUTIONS
108; ==========
109;
110; - fplod 20110809T110911Z aedon.locean-ipsl.upmc.fr (Darwin)
111;
112;   * complete descritption
113;   * remove v* from filenames (in and out)
114;   * usage of ${PROJECT_OD}
115;   * remove return statement
116;   * add test on IO files
117;
118; - fplod 20101217T140745Z aedon.locean-ipsl.upmc.fr (Darwin)
119;
120;   * remove hard coded directory for mask_oaflux_30N30S.nc
121;
122; - fplod 20101214T112131Z aedon.locean-ipsl.upmc.fr (Darwin)
123;
124;   * add graph
125;
126; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin)
127;
128;   * minimal header
129;
130; - pbk 2008
131;
132;   * creation
133;
134;-
135;
136pro TropFlux_19890101_20091231
137;
138@cm_4cal
139@cm_4data
140@cm_4mesh
141@cm_4data
142@cm_project
143;
144; check for input directory
145;
146; test if ${PROJECT_ID} defined
147CASE project_id_env OF
148    ''  :  BEGIN
149     msg = 'eee : ${PROJECT_ID} is not defined'
150     ras = report(msg)
151     STOP
152           END
153 ELSE: BEGIN
154     msg = 'iii : ${PROJECT_ID} is ' + project_id_env
155     ras = report(msg)
156       END
157ENDCASE
158;
159iodirin = isadirectory(project_id_env)
160;
161; existence and protection of ${PROJECT_ID}
162IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
163   msg = 'eee : the directory' + iodirin  + ' is not accessible.'
164   ras = report(msg)
165   STOP
166ENDIF
167;
168; build mask filename
169filename_msk='mask_oaflux_30N30S.nc'
170;
171; check if this file exists
172msg='iii : looking for ' + filename_msk
173ras = report(msg)
174fullfilename_msk = isafile(iodirin + filename_msk, NEW=0, /MUST_EXIST)
175IF fullfilename_msk[0] EQ '' THEN BEGIN
176   msg = 'eee : the file ' + fullfilename_msk + ' was not found.'
177   ras = report(msg)
178   STOP
179ENDIF
180;
181; test if ${PROJECT_OD} defined
182CASE project_od_env OF
183  '' : BEGIN
184         msg = 'eee : ${PROJECT_OD} is not defined'
185         ras = report(msg)
186       STOP
187       END
188  ELSE: BEGIN
189          msg = 'iii : ${PROJECT_OD} is ' + project_od_env
190          ras = report(msg)
191        END
192 ENDCASE
193;
194; check if output data will be possible
195iodirout = isadirectory(project_od_env)
196;
197; existence and protection for reading
198IF (FILE_TEST(iodirout, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
199   msg = 'eee : the directory' + iodirout  + ' is not accessible.'
200   ras = report(msg)
201   STOP
202ENDIF
203;
204; existence and protection for writing
205IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN
206    msg = 'eee : the directory' + iodirout  + ' was not found.'
207    ras = report(msg)
208    STOP
209ENDIF
210;
211; build sst data filename
212filename_sst='TropFlux_sst_19890101_20091231.nc'
213;
214; check if this file exists
215msg='iii : looking for ' + filename_sst
216ras = report(msg)
217fullfilename_sst = isafile(iodirout + filename_sst, NEW=0, /MUST_EXIST)
218IF fullfilename_sst[0] EQ '' THEN BEGIN
219   msg = 'eee : the file ' + fullfilename_sst + ' was not found.'
220   ras = report(msg)
221   STOP
222ENDIF
223;
224; build ws data filename
225filename_ws='TropFlux_ws_19890101_20091231.nc'
226;
227; check if this file exists
228msg='iii : looking for ' + filename_ws
229ras = report(msg)
230fullfilename_ws = isafile(iodirout + filename_ws, NEW=0, /MUST_EXIST)
231IF fullfilename_ws[0] EQ '' THEN BEGIN
232   msg = 'eee : the file ' + fullfilename_ws + ' was not found.'
233   ras = report(msg)
234   STOP
235ENDIF
236;
237; build swr data filename
238filename_swr='TropFlux_swr_19890101_20091231_BLND.nc'
239;
240; check if this file exists
241msg='iii : looking for ' + filename_swr
242ras = report(msg)
243fullfilename_swr = isafile(iodirout + filename_swr, NEW=0, /MUST_EXIST)
244IF fullfilename_swr[0] EQ '' THEN BEGIN
245   msg = 'eee : the file ' + fullfilename_swr + ' was not found.'
246   ras = report(msg)
247   STOP
248ENDIF
249;
250; build lwr data filename
251filename_lwr='TropFlux_lwr_19890101_20091231.nc'
252;
253; check if this file exists
254msg='iii : looking for ' + filename_lwr
255ras = report(msg)
256fullfilename_lwr = isafile(iodirout + filename_lwr, NEW=0, /MUST_EXIST)
257IF fullfilename_lwr[0] EQ '' THEN BEGIN
258   msg = 'eee : the file ' + fullfilename_lwr + ' was not found.'
259   ras = report(msg)
260   STOP
261ENDIF
262;
263; build t2m data filename
264filename_t2m='TropFlux_t2m_19890101_20091231.nc'
265;
266; check if this file exists
267msg='iii : looking for ' + filename_t2m
268ras = report(msg)
269fullfilename_t2m = isafile(iodirout + filename_t2m, NEW=0, /MUST_EXIST)
270IF fullfilename_t2m[0] EQ '' THEN BEGIN
271   msg = 'eee : the file ' + fullfilename_t2m + ' was not found.'
272   ras = report(msg)
273   STOP
274ENDIF
275;
276; build q2m data filename
277filename_q2m='TropFlux_q2m_19890101_20091231.nc'
278;
279; check if this file exists
280msg='iii : looking for ' + filename_q2m
281ras = report(msg)
282fullfilename_q2m = isafile(iodirout + filename_q2m, NEW=0, /MUST_EXIST)
283IF fullfilename_q2m[0] EQ '' THEN BEGIN
284   msg = 'eee : the file ' + fullfilename_q2m + ' was not found.'
285   ras = report(msg)
286   STOP
287ENDIF
288;
289; build wg data filename
290filename_wg='TropFlux_gustiness_19890101_20091231.nc'
291;
292; check if this file exists
293msg='iii : looking for ' + filename_wg
294ras = report(msg)
295fullfilename_wg = isafile(iodirout + filename_wg, NEW=0, /MUST_EXIST)
296IF fullfilename_wg[0] EQ '' THEN BEGIN
297   msg = 'eee : the file ' + fullfilename_wg + ' was not found.'
298   ras = report(msg)
299   STOP
300ENDIF
301;
302da1=19880101 & da2=20101231
303;
304initncdf, fullfilename_msk
305msk=ncdf_lec(fullfilename_msk,var='msk')
306
307initncdf, fullfilename_sst
308
309ws=read_ncdf('ws',da1,da2,file=fullfilename_ws,/nostr)
310wg=read_ncdf('wg',da1-1,da2,file=fullfilename_wg,/nostr)
311tt=time & jpt=n_elements(time)
312sst=read_ncdf('sst',da1,da2,file=fullfilename_sst,/nostr)
313sst=reform(sst-273.15)
314swd=read_ncdf('swr',da1,da2,file=fullfilename_swr,/nostr)
315lw=read_ncdf('lwr',da1,da2,file=fullfilename_lwr,/nostr)
316swd=swd/0.94  ;; converting from net swr to downward swr
317
318t2m=read_ncdf('t2m',da1,da2,file=fullfilename_t2m,/nostr)-273.15   ; in C
319q2m=read_ncdf('q2m',da1,da2,file=fullfilename_q2m,/nostr)    ; in g/kg
320
321w=sqrt(ws*ws+wg*wg)  ;; wind corrected for gustiness
322w=ws
323tmask=msk
324help, ws,w,wg,u,sst,swd,t2m,q2m
325ocean=where(msk eq 1,compl=land)
326valmask=1.e20
327time=tt & jpt=n_elements(time)
328
329;
330;; Constants for flux computation
331;
332zu=10.                          ; height of wind speed measurement (m)
333us=0.                           ; surf current (m/s)
334zt=2.                           ; Height of air T measurement (m)
335zq=2.                           ; height of humidity measurement (m)
336P=1008.                         ; Pressure
337zi=600.                         ; Inversion height (m)
338jcool=0                         ; Compute cool-skin
339jwave=0                         ; No waves
340twave=5.
341hwave=1.
342
343caldat, time,mon,day,yea
344swr=fltarr(jpi,jpj,jpt)+1.e20
345lwr=fltarr(jpi,jpj,jpt)+1.e20
346lat=fltarr(jpi,jpj,jpt)+1.e20
347sen=fltarr(jpi,jpj,jpt)+1.e20
348lwnet_clrk=fltarr(jpi,jpj,jpt)+1.e20
349;Ch=fltarr(jpi,jpj,jpt)+1.e20
350;Ce=fltarr(jpi,jpj,jpt)+1.e20
351junk=fltarr(jpi,jpj,jpt)+1.e20
352
353for jt=0,jpt-1 do begin
354  jday=time(jt)-julday(1,1,yea(jt))
355  print, 'Computing Fluxes ',jt,' / ',jpt-1
356;
357;  P=msl(*,*,jt) & P=P(ocean)
358  wn=w(*,*,jt) & wn=wn(ocean)          ; wind speed (m/s)
359  ts=sst(*,*,jt) & ts=ts(ocean)        ; Bulk sst (°C)
360  t=t2m(*,*,jt) & t=t(ocean)           ; 2m Air T (°C)
361  qs=qsee(ts,P)                        ; Sea surface sat. spec. humidity (g/kg)
362  q=q2m(*,*,jt) & q=q(ocean)           ; 2m AIr specific humidity  (g/kg)
363  Rs=swd(*,*,jt) & Rs=Rs(ocean)        ; Downward solar flux (W/m2)
364  ylat=gphit(ocean)
365;  cld=calc_cloud(jday,Rs,ylat)
366  cld=calc_cloud_vlat(jday,Rs,ylat)
367;  Rl=lwdown_clark(ts,q,cld,t,P)        ; Downward IR flux (W/m2)
368  Rl=lw(*,*,jt) & Rl=Rl(ocean)
369  rain=0.
370  lw_clrk=-lwnet_clark(ts,q,cld,t,P)
371;  junk(*,*,jt)=lw_clrk
372;
373;stop
374  y=cor30a(wn,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,ylat,jcool,jwave,twave,hwave)
375;
376
377; A few punctual missing values (coare does not converge): filled by spatial extrapolation
378  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y(*,0)) & tab(ocean)=x & m=finite(tab) & tab=extrapolate(tab,m) & swr(*,*,jt)=tab*msk+valmask*(1-msk)
379  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y(*,1)) & tab(ocean)=x & m=finite(tab) & tab=extrapolate(tab,m) & lwr(*,*,jt)=tab*msk+valmask*(1-msk)
380  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y(*,2)) & tab(ocean)=x & m=finite(tab) & tab=extrapolate(tab,m) & lat(*,*,jt)=tab*msk+valmask*(1-msk)
381  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y(*,3)) & tab(ocean)=x & m=finite(tab) & tab=extrapolate(tab,m) & sen(*,*,jt)=tab*msk+valmask*(1-msk)
382;  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(lw_clrk) & tab(ocean)=x & m=finite(tab) & tab=extrapolate(tab,m) & lwnet_clrk(*,*,jt)=tab*msk+valmask*(1-msk)
383;  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y(*,6)) & tab(ocean)=x & m=finite(tab) & tab=extrapolate(tab,m) & Ch(*,*,jt)=tab*msk+valmask*(1-msk)
384;  tab=fltarr(jpi,jpj)+!values.f_nan & x=reform(y(*,7)) & tab(ocean)=x & m=finite(tab) & tab=extrapolate(tab,m) & Ce(*,*,jt)=tab*msk+valmask*(1-msk)
385endfor
386tt=time
387time=timegen(7670, start=julday(1,1,1989,0), units='days') & jpt=n_elements(time)
388;
389cda0=string(jul2date(time(0)),format='(i8.8)')
390cda1=string(jul2date(time(jpt-1)),format='(i8.8)')
391tt=time-julday(1,1,1950,00,00,00)
392xlon=reform(glamt(*,0) ) & ylat=reform(gphit(0,*))
393
394ncfile='!${PROJECT_OD}/TropFlux_19890101_20091231.nc'
395lon_attr={units:'degrees_east',long_name:'Longitude'}
396lat_attr={units:'degrees_north',long_name:'Latitude'}
397swr_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net shortwave radiation',short_name:'swr',axis:'TYX'}
398lwr_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net longwave radiation',short_name:'lwr',axis:'TYX'}
399lwr_clrk_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net longwave radiation from clark',short_name:'lwr',axis:'TYX'}
400lhf_attr={units:'W/m2',missing_value:valmask,long_name:'Surface latent flux',short_name:'lhf',axis:'TYX'}
401shf_attr={units:'W/m2',missing_value:valmask,long_name:'Surface sensible flux',short_name:'shf',axis:'TYX'}
402time_attr={units:'days since 1950-01-01 00:00:00',long_name:'Time axis',time_origin:' 1950-JAN-01 00:00:00'}
403Ch_attr={units:'',missing_value:valmask,long_name:'heat transfer coefficient at zt',short_name:'Ch',axis:'TYX'}
404Ce_attr={units:'',missing_value:valmask,long_name:'moisture transfer coefficient at zq',short_name:'Ce',axis:'TYX'}
405
406globattr={source:'Fluxes for the Global Tropical Ocean - TropFlux',timerange:cda0+' - '+cda1}
407
408help, swr,lwr,lat,sen,tt,xlon,ylat
409
410ncfields = 'swr[longitude,latitude,time]=swr:swr_attr; ' $
411;          +'lwr_coare[longitude,latitude,time]=lwr:lwr_attr; ' $
412          +'lwr[longitude,latitude,time]=lwnet_clrk:lwr_clrk_attr; ' $
413          +'lhf[longitude,latitude,time]=lat:lhf_attr; ' $
414          +'shf[longitude,latitude,time]=sen:shf_attr; ' $
415;          +'Ch[longitude,latitude,time]=Ch:Ch_attr; ' $
416;          +'Ce[longitude,latitude,time]=Ce:Ce_attr; ' $
417                      + 'longitude[]=xlon:lon_attr; ' $
418                      + 'latitude[]=ylat:lat_attr; ' $
419                      + 'tt[*time]=tt:time_attr ' $
420                      + ' @ globattr'
421
422@ncdf_quickwrite
423
424end
Note: See TracBrowser for help on using the repository browser.