source: trunk/src/TropFlux_19890101_20091231.pro @ 70

Last change on this file since 70 was 50, checked in by pinsard, 13 years ago

replace TROPFLUX by PROJECT
  • Property svn:executable set to *
File size: 10.1 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;
17;     .. graphviz::
18;
19;        digraph TropFlux_19890101_20091231 {
20;           graph [
21;           rankdir="LR",
22;           ]
23;           mask [shape=ellipse,fontname=Courier,label="${PROJECT_ID}/mask_oaflux_30N30S.nc"];
24;           file_sst [shape=ellipse,fontname=Courier,label="/Volumes/Iomega_HDD/TropFlux/input_cor/full_cor/TropFlux_sst_19890101_20091231_v20.nc"];
25;           file_wind [shape=ellipse,fontname=Courier,label="/Volumes/Iomega_HDD/TropFlux/input_cor/full_cor/TropFlux_ws_19890101_20091231_v20.nc"];
26;           file_wg [shape=ellipse,fontname=Courier,label="/Users/pkb/data/TropFlux/TropFlux_gustiness_19890101_20091231_v50.nc"];
27;           file_sw [shape=ellipse,fontname=Courier,label="/Users/pkb/data/TropFlux/TropFlux_swr_19890101_20091231_BLND_v50.nc"];
28;           file_lw [shape=ellipse,fontname=Courier,label="/Volumes/Iomega_HDD/TropFlux/input_cor/full_cor/TropFlux_lwr_19890101_20091231_v2.nc"];
29;           file_air [shape=ellipse,fontname=Courier,label="/Users/pkb/data/TropFlux/TropFlux_t2m_19890101_20091231_v50.nc"];
30;           file_q [shape=ellipse,fontname=Courier,label="/Volumes/Iomega_HDD/TropFlux/input_cor/full_cor/TropFlux_q2m_19890101_20091231_v20.nc"];
31;
32;           ncfile [shape=ellipse,fontname=Courier,label="/Users/pkb/data/TropFlux/TropFlux_19890101_20091231_v51.nc"];
33;
34;           TropFlux_19890101_20091231 [shape=box,
35;           fontname=Courier,
36;           color=blue,
37;           URL="http://forge.ipsl.jussieu.fr/tropflux/browser/trunk/src/TropFlux_19890101_20091231.pro",
38;           label="${PROJECT}/src/TropFlux_19890101_20091231.pro"];
39;
40;           {mask file_sst file_wind file_wg file_sw file_lw file_air file_q} -> {TropFlux_19890101_20091231} -> {ncfile}
41;
42;          }
43;
44; SEE ALSO
45; ========
46;
47; :ref:`project_profile.sh`
48;
49; :func:`report <saxo:report>`
50; :func:`isadirectory <saxo:isadirectory>`
51; :func:`isafile <saxo:isafile>`
52; :func:`initncdf <saxo:initncdf>`
53; :func:`ncdf_lec <saxo:ncdf_lec>`
54; :func:`read_ncdf <saxo:read_ncdf>`
55; :func:`caldat <saxo:caldat>`
56; :func:`julday <saxo:julday>`
57; :func:`ncdf_quickwrite <saxo:ncdf_quickwrite>`
58;
59; :func:`cor30a`
60;
61; EXAMPLES
62; ========
63;
64; ::
65;
66;  IDL> tropflux_19890101_20091231
67;
68; TODO
69; ====
70;
71; hard coded directory - usage of ${PROJECT_ID}
72;
73; coding rules
74;
75; KNOWN ISSUES
76; ============
77;
78; test of existence of fullfilename_msk not very efficient because
79; MUST_EXIST keyword of :func:`isafile <saxo:isafile>` not yet implemented
80;
81; EVOLUTIONS
82; ==========
83;
84; - fplod 20101217T140745Z aedon.locean-ipsl.upmc.fr (Darwin)
85;
86;   * remove hard coded directory for mask_oaflux_30N30S.nc
87;
88; - fplod 20101214T112131Z aedon.locean-ipsl.upmc.fr (Darwin)
89;
90;   * add graph
91;
92; - fplod 20101214T093615Z aedon.locean-ipsl.upmc.fr (Darwin)
93;
94;   * minimal header
95;
96; - pbk 2008
97;
98;   * creation
99;
100;-
101;
102pro TropFlux_19890101_20091231
103@common
104@cm_project
105;
106; check for input directory
107;
108; test if ${PROJECT_ID} defined
109CASE project_id_env OF
110    ''  :  BEGIN
111     msg = 'eee : ${PROJECT_ID} is not defined'
112     ras = report(msg)
113     STOP
114           END
115 ELSE: BEGIN
116     msg = 'iii : ${PROJECT_ID} is ' + project_id_env
117     ras = report(msg)
118       END
119ENDCASE
120;
121iodirin = isadirectory(project_id_env)
122;
123; existence and protection of ${PROJECT_ID}
124IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
125   msg = 'eee : the directory' + iodirin  + ' is not accessible.'
126   ras = report(msg)
127   STOP
128ENDIF
129;
130; build mask filename
131filename_msk='mask_oaflux_30N30S.nc'
132;
133; check if this file exists
134fullfilename_msk = isafile(iodirin + filename_msk, NEW=0, /MUST_EXIST)
135IF fullfilename_msk[0] EQ '' THEN BEGIN
136   msg = 'eee : the file ' + fullfilename_msk + ' was not found.'
137   ras = report(msg)
138   STOP
139ENDIF
140
141da1=19880101 & da2=20101231
142;
143initncdf, fullfilename_msk
144msk=ncdf_lec(fullfilename_msk,var='msk')
145
146dir='/Volumes/Iomega_HDD/TropFlux/input_cor/full_cor/'
147dir1='/Users/pkb/data/TropFlux/'
148
149file_sst=dir+'TropFlux_sst_19890101_20091231_v20.nc'
150file_wind=dir+'TropFlux_ws_19890101_20091231_v20.nc'
151file_sw=dir1+'TropFlux_swr_19890101_20091231_BLND_v50.nc'
152file_lw=dir+'TropFlux_lwr_19890101_20091231_v2.nc'
153file_air=dir1+'TropFlux_t2m_19890101_20091231_v50.nc'
154file_q=dir+'TropFlux_q2m_19890101_20091231_v20.nc'
155file_wg=dir1+'TropFlux_gustiness_19890101_20091231_v50.nc'
156
157initncdf, file_sst
158
159ws=read_ncdf('ws',da1,da2,file=file_wind,/nostr)
160wg=read_ncdf('wg',da1-1,da2,file=file_wg,/nostr)
161tt=time & jpt=n_elements(time)
162sst=read_ncdf('sst',da1,da2,file=file_sst,/nostr)
163sst=reform(sst-273.15)
164swd=read_ncdf('swr',da1,da2,file=file_sw,/nostr)
165lw=read_ncdf('lwr',da1,da2,file=file_lw,/nostr)
166swd=swd/0.94  ;; converting from net swr to downward swr
167
168t2m=read_ncdf('t2m',da1,da2,file=file_air,/nostr)-273.15   ; in C
169q2m=read_ncdf('q2m',da1,da2,file=file_q,/nostr)    ; in g/kg
170
171w=sqrt(ws*ws+wg*wg)  ;; wind corrected for gustiness
172w=ws
173tmask=msk
174help, ws,w,wg,u,sst,swd,t2m,q2m
175ocean=where(msk eq 1,compl=land)
176valmask=1.e20
177time=tt & jpt=n_elements(time)
178
179;
180;; Constants for flux computation
181;
182zu=10.                          ; height of wind speed measurement (m)
183us=0.                           ; surf current (m/s)
184zt=2.                           ; Height of air T measurement (m)
185zq=2.                           ; height of humidity measurement (m)
186P=1008.                         ; Pressure
187zi=600.                         ; Inversion height (m)
188jcool=0                         ; Compute cool-skin
189jwave=0                         ; No waves
190twave=5.
191hwave=1.
192
193caldat, time,mon,day,yea
194swr=fltarr(jpi,jpj,jpt)+1.e20
195lwr=fltarr(jpi,jpj,jpt)+1.e20
196lat=fltarr(jpi,jpj,jpt)+1.e20
197sen=fltarr(jpi,jpj,jpt)+1.e20
198lwnet_clrk=fltarr(jpi,jpj,jpt)+1.e20
199;Ch=fltarr(jpi,jpj,jpt)+1.e20
200;Ce=fltarr(jpi,jpj,jpt)+1.e20
201junk=fltarr(jpi,jpj,jpt)+1.e20
202
203for jt=0,jpt-1 do begin
204  jday=time(jt)-julday(1,1,yea(jt))
205  print, 'Computing Fluxes ',jt,' / ',jpt-1
206;
207;  P=msl(*,*,jt) & P=P(ocean)
208  wn=w(*,*,jt) & wn=wn(ocean)          ; wind speed (m/s)
209  ts=sst(*,*,jt) & ts=ts(ocean)        ; Bulk sst (°C)
210  t=t2m(*,*,jt) & t=t(ocean)           ; 2m Air T (°C)
211  qs=qsee(ts,P)                        ; Sea surface sat. spec. humidity (g/kg)
212  q=q2m(*,*,jt) & q=q(ocean)           ; 2m AIr specific humidity  (g/kg)
213  Rs=swd(*,*,jt) & Rs=Rs(ocean)        ; Downward solar flux (W/m2)
214  ylat=gphit(ocean)
215;  cld=calc_cloud(jday,Rs,ylat)
216  cld=calc_cloud_vlat(jday,Rs,ylat)
217;  Rl=lwdown_clark(ts,q,cld,t,P)        ; Downward IR flux (W/m2)
218  Rl=lw(*,*,jt) & Rl=Rl(ocean)
219  rain=0.
220  lw_clrk=-lwnet_clark(ts,q,cld,t,P)
221;  junk(*,*,jt)=lw_clrk
222;
223;stop
224  y=cor30a(wn,us,ts,t,Qs,Q,Rs,Rl,rain,zi,P,zu,zt,zq,ylat,jcool,jwave,twave,hwave)
225;
226
227; A few punctual missing values (coare does not converge): filled by spatial extrapolation
228  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)
229  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)
230  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)
231  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)
232;  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)
233;  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)
234;  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)
235endfor
236tt=time
237time=timegen(7670, start=julday(1,1,1989,0), units='days') & jpt=n_elements(time)
238;
239cda0=string(jul2date(time(0)),format='(i8.8)')
240cda1=string(jul2date(time(jpt-1)),format='(i8.8)')
241tt=time-julday(1,1,1950,00,00,00)
242xlon=reform(glamt(*,0) ) & ylat=reform(gphit(0,*))
243
244ncfile='!/Users/pkb/data/TropFlux/TropFlux_19890101_20091231_v51.nc'
245lon_attr={units:'degrees_east',long_name:'Longitude'}
246lat_attr={units:'degrees_north',long_name:'Latitude'}
247swr_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net shortwave radiation',short_name:'swr',axis:'TYX'}
248lwr_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net longwave radiation',short_name:'lwr',axis:'TYX'}
249lwr_clrk_attr={units:'W/m2',missing_value:valmask,long_name:'Surface net longwave radiation from clark',short_name:'lwr',axis:'TYX'}
250lhf_attr={units:'W/m2',missing_value:valmask,long_name:'Surface latent flux',short_name:'lhf',axis:'TYX'}
251shf_attr={units:'W/m2',missing_value:valmask,long_name:'Surface sensible flux',short_name:'shf',axis:'TYX'}
252time_attr={units:'days since 1950-01-01 00:00:00',long_name:'Time axis',time_origin:' 1950-JAN-01 00:00:00'}
253Ch_attr={units:'',missing_value:valmask,long_name:'heat transfer coefficient at zt',short_name:'Ch',axis:'TYX'}
254Ce_attr={units:'',missing_value:valmask,long_name:'moisture transfer coefficient at zq',short_name:'Ce',axis:'TYX'}
255
256globattr={source:'Fluxes for the Global Tropical Ocean - TropFlux',timerange:cda0+' - '+cda1}
257
258help, swr,lwr,lat,sen,tt,xlon,ylat
259
260ncfields = 'swr[longitude,latitude,time]=swr:swr_attr; ' $
261;          +'lwr_coare[longitude,latitude,time]=lwr:lwr_attr; ' $
262          +'lwr[longitude,latitude,time]=lwnet_clrk:lwr_clrk_attr; ' $
263          +'lhf[longitude,latitude,time]=lat:lhf_attr; ' $
264          +'shf[longitude,latitude,time]=sen:shf_attr; ' $
265;          +'Ch[longitude,latitude,time]=Ch:Ch_attr; ' $
266;          +'Ce[longitude,latitude,time]=Ce:Ce_attr; ' $
267                      + 'longitude[]=xlon:lon_attr; ' $
268                      + 'latitude[]=ylat:lat_attr; ' $
269                      + 'tt[*time]=tt:time_attr ' $
270                      + ' @ globattr'
271
272@ncdf_quickwrite
273
274return
275end
Note: See TracBrowser for help on using the repository browser.