Hi, all. I am trying to calculate longwave radiation with ground-
measured air temperature and relative humidity data, but my output
file looks crazy. Is there anyone can help me?
What I've done so far is as follows.
Please take a look, and leave any comments. Thanks.
-------------------------------
Pro CALC_Longwave_NWS
; ********************************** Get data from NWS
**********************************************************
s_time=systime(1)
;year=2002
WorkSp = 'D:\MODIS07\Longwave_Workspace\'
WorkSpOut = 'D:\MODIS07\Longwave_Workspace\output\'
; sitefile = strcompress(WorkSp+'site_lon_lat.txt', /remove_all)
;site_id year month day h01 h02 h03 h04 h05 h06 h07 h08 h09 h10 h11
h12 h13 h14 h15 h16 h17 h18 h19 h20 h21 h22 h23 h24
;90 2002 1 1 25 25 26 26 26 25 26
28 31 34 36 39 40 38 38 37 35 33 35 37 34
34 35 29
;90 2002 1 2 30 32 29 28 31 33 55
65 74 68 72 75 71 83 85 77 69 60 51 43 41
34 34 29
NWS_temp = strcompress(WorkSp+'temp2002.txt', /remove_all) ;
data1(28, nrow), temperature in 2002
nrow1 = file_lines(NWS_temp)
data1 = fltarr(28, nrow1)
close,/all
Openr, 1, NWS_temp
readf, 1, data1
close, 1
cd, 'D:\MODIS07\Longwave_Workspace\'
readcol, 'temp2002.txt', site_id, year, month, day, h01, h02, h03,
h04, h05, h06, h07, h08, h09, h10, h11, h12, h13, h14, h15, h16, h17,
h18, h19, h20, h21, h22, h23, h24,
Format='(I,I,I,I,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F ,F,F)'
NWS_hum = strcompress(WorkSp+'hum2002.txt' , /remove_all) ; data2(278,
nrow-1), humidity in 2002
;nrow2 = file_lines(NWS_hum)
data2 = fltarr(28, nrow1)
close,/all
Openr, 2, NWS_hum
readf, 2, data2
close, 2
readcol, 'temp2002.txt', site_idh, yearh, monthh, dayh, h01h, h02h,
h03h, h04h, h05h, h06h, h07h, h08h, h09h, h10h, h11h, h12h, h13h,
h14h, h15h, h16h, h17h, h18h, h19h, h20h, h21h, h22h, h23h, h24h,
Format='(I,I,I,I,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F ,F,F)'
Es = fltarr(28, nrow1) ; Saturated vapor pressure, from NWS,
Clausius-Clapeyron equation
Ea = fltarr(28, nrow1) ; actual vapor pressure, from NWS
Brunt_temp = fltarr(28, nrow1) ; DownRl, Brunt, 1932
Brunt = fltarr(28, nrow1) ; DownRl, Brunt, 1932
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;read input data for each site
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
; calculating Ea,actual vapor pressure
Lv = 2.45 * 10^6 ; the latent heat of vaporization, J /kg
Rv = 461.0 ; the gas constant for water vapor, J/ kg/K
sigma = 5.67*10.0^(-8) ; Steffan-Boltzmann constant, W/m^2/K^(-4)
Es = 611*exp((Lv/Rv)*((1/273.16)-(1/data1)))
Ea = data2*Es/100
; ************** Downward Longwave Radiation 2-Brunt, 1932
***************************
a1 = 0.605
b1 = 0.048
Brunt_temp = (a1 + b1*Ea^0.5)*Sigma*data1^4
;Brunt_out = strcompress(WorkSpOut+'Bruntout2002.txt', /remove_all)
bsite_id=fltarr(nrow1)
byear=fltarr(nrow1)
bmonth=fltarr(nrow1)
bday=fltarr(nrow1)
bh01=fltarr(nrow1)
bh02=fltarr(nrow1)
bh03=fltarr(nrow1)
bh04=fltarr(nrow1)
bh05=fltarr(nrow1)
bh06=fltarr(nrow1)
bh07=fltarr(nrow1)
bh08=fltarr(nrow1)
bh09=fltarr(nrow1)
bh10=fltarr(nrow1)
bh11=fltarr(nrow1)
bh12=fltarr(nrow1)
bh13=fltarr(nrow1)
bh14=fltarr(nrow1)
bh15=fltarr(nrow1)
bh16=fltarr(nrow1)
bh17=fltarr(nrow1)
bh18=fltarr(nrow1)
bh19=fltarr(nrow1)
bh20=fltarr(nrow1)
bh21=fltarr(nrow1)
bh22=fltarr(nrow1)
bh23=fltarr(nrow1)
bh24=fltarr(nrow1)
openw, lun, 'brunt_temp.txt', /get_lun
printf, lun, Brunt_temp, bsite_id, byear, bmonth, bday, bh01, bh02,
bh03, bh04, bh05, bh06, bh07, bh08, bh09, bh10, bh11, bh12, bh13,
bh14, bh15, bh16, bh17, bh18, bh19, $
bh20, bh21, bh22, bh23, bh24, format='(4I7, x, 24(F7.2))'
close, lun
;readcol, 'brunt_temp.txt', bsite, byear, bmonth, bday, b01, b02, b03,
b04, b05, b06, b07, b08, b09, b10, b11, b12, b13, b14, b15, b16, b17,
b18, b19, b20, b21, b22, b23, b24, $
;
Format='(I,I,I,I,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F,F ,F,F)'
openw, lun2, 'brunt_temp.txt', /get_lun
;for i=0, nrow1-1 do begin
for i=0, 9 do begin
print, 'Now I am writing ', i, ' th line !!!'
printf, lun2, site_id, year, month, day, bh01, bh02, bh03, bh04,
bh05, bh06, bh07, bh08, bh09, bh10, bh11, bh12, bh13, bh14, bh15,
bh16, bh17, bh18, bh19, bh20, bh21, bh22, bh23, bh24, format='(4I7, x,
24(F7.2))'
endfor
close, lun2
f_time = systime(1)
etime = f_time-s_time
etime_str = string(etime, format='(F7.2)')
print, 'Elapsed time for this process = ', etime_str
end
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
