; Filename: read_qscat3.pro ; ; Usage: ; ; To run this programs, use the following command: ; ; IDL> read_qscat3 ; ; Description: ; ; This file contains one (1) IDL procedure and three (3) functions ; used to read the QuikSCAT Level 3 data in Hierarchical Data ; Format (HDF). The routines are as follows. ; ; 1. correct_uint16: a function to correct variables ; stored as 16-bit unsigned integers. ; ; 2. correct_uint32: a function to correct variables ; stored as 32-bit unsigned integers. ; ; 3. get_sds: a function to read the data which are ; stored as Scientific Data Sets. ; ; 4. read_qscat3: the main procedure. Calls get_vdata and ; get_sds. ; ; Notes: ; ; 1. The directory on your local system which contains the QuikSCAT ; L3 data must be input in program read_qscat3.pro. (Search on ; "data_dir" to find this line quickly.) ; ; 2. The L3 data are read in their entirety. If you would like ; to read these data in slabs, change the value of "slab_size". ; ; 3. The correction of Unsigned Integers is currently hard-coded for ; versions of IDL prior to 5.2. ; ; 4. Please send all comments and questions concerning these routines ; to qscat@podaac.jpl.nasa.gov. ; ; ; 10/4/99 K.L. Perry, C.S. Hsu, R.S. Dunbar ; 12/10/03 T. Lungu ; Revisions: ; ; 5/5/2000 Updated code to read asc_time_frac, des_time_frac, ; asc_rain_prob, des_rain_prob, asc_rain_flag, and ; des_rain_flag. K.L. Perry ; 12/10/03 Replaced pickfile by dialog_pickfile, since the former is ; no longer supported by IDL 6.0. ;================================================================== ; Copyright (c) 1999, California Institute of Technology ;================================================================== ;================================================================== ; correct_uint16: a function to correct data which are stored ; as 16-bit unsigned integers for versions of ; IDL prior to 5.2 ;================================================================== function correct_uint16,sds_data w=where(sds_data lt 0,cnt) if (cnt gt 0) then begin sds_data(w)=sds_data(w)+65536. w2=where(sds_data lt 0,cnt2) endif return,sds_data end ;================================================================== ; correct_uint32: a function to correct data which are stored ; as 32-bit unsigned integers for versions of ; IDL prior to 5.2 ;================================================================== function correct_uint32,sds_data w=where(sds_data lt 0,cnt) if (cnt gt 0) then begin sds_data(w)=sds_data(w)+4294967296. w2=where(sds_data lt 0,cnt2) endif return,sds_data end ;================================================================== ; get_sds: a function to read the data which are stored ; as Scientific Data Sets. This function also ; multiplies the data by the calibration and ; subtracts the offset. ;================================================================== function get_sds,sd_id,sds_name,slab_start,slab_size index=HDF_SD_NAMETOINDEX(sd_id,sds_name) sds_id=HDF_SD_SELECT(sd_id,index) HDF_SD_GETINFO,sds_id,ndims=ndims,dims=dims,caldata=cal,type=data_type ; make sure edge, start and stride have correct array sizes edge=intarr(ndims) start=intarr(ndims) stride=intarr(ndims) for i=0,ndims-1 do begin edge(i)=dims(i) start(i)=0 stride(i)=1 endfor edge(ndims-1)=slab_size start(ndims-1)=slab_start HDF_SD_GETDATA,sds_id,data,stride=stride,start=start,count=edge ;Correct Unsigned Integers ;; note: Versions of IDL prior to 5.2 do not handle unsigned ;; integers properly. These versions also do not identify ;; DFNT_UINT's using HDF_SD_GETINFO. Therefore, the ;; following hard code will be included until IDL 5.2 is ;; "standard". --KLP. ;; ;; UINT16 if ((sds_name eq "asc_avg_wind_speed") or $ (sds_name eq "des_avg_wind_speed") or $ (sds_name eq "asc_time_frac") or $ (sds_name eq "des_time_frac") or $ (sds_name eq "asc_rain_prob") or $ (sds_name eq "des_rain_prob")) then begin data=correct_uint16(float(data)) endif ;; UINT32 if ((sds_name eq "asc_avg_wind_speed_sq") or $ (sds_name eq "des_avg_wind_speed_sq")) then begin data=correct_uint16(float(data)) endif ; Apply the scale and offset rdata = data*cal.cal - cal.offset HDF_SD_ENDACCESS,sds_id return,rdata end ;================================================================== ; read_qscat3: the main procedure in this file. Calls are made ; to get_sds to read the data. The results are then ; printed to the screen. ;================================================================== pro read_qscat3 ;Select a Level 3 file ;***** Change data_dir to suit your local system data_dir="./" filename=dialog_pickfile(filter='QS_XWGRD3*', PATH = data_dir, /READ) print,' ' print,'FILENAME: ',filename print,' ' ;Read the Scientific Data Sets sd_id=HDF_SD_START(filename,/READ) ;Get the x-dimension in order to read an entire SDS index=HDF_SD_NAMETOINDEX(sd_id,'asc_wvc_count') sds_id=HDF_SD_SELECT(sd_id,index) HDF_SD_GETINFO,sds_id,ndims=rank,dims=dims,label=name, $ type=data_type,caldata=cal slab_size=strtrim(dims(1),2) ir=0 ; Read the SDSs asc_avg_wind_speed= get_sds(sd_id,'asc_avg_wind_speed',ir,slab_size) des_avg_wind_speed= get_sds(sd_id,'des_avg_wind_speed',ir,slab_size) asc_avg_wind_vel_u= get_sds(sd_id,'asc_avg_wind_vel_u',ir,slab_size) des_avg_wind_vel_u= get_sds(sd_id,'des_avg_wind_vel_u',ir,slab_size) asc_avg_wind_vel_v= get_sds(sd_id,'asc_avg_wind_vel_v',ir,slab_size) des_avg_wind_vel_v= get_sds(sd_id,'des_avg_wind_vel_v',ir,slab_size) asc_avg_wind_speed_sq= get_sds(sd_id,'asc_avg_wind_speed_sq',ir,slab_size) des_avg_wind_speed_sq= get_sds(sd_id,'des_avg_wind_speed_sq',ir,slab_size) asc_wvc_count= get_sds(sd_id,'asc_wvc_count',ir,slab_size) des_wvc_count= get_sds(sd_id,'des_wvc_count',ir,slab_size) asc_time_frac= get_sds(sd_id,'asc_time_frac',ir,slab_size) des_time_frac= get_sds(sd_id,'des_time_frac',ir,slab_size) asc_rain_prob= get_sds(sd_id,'asc_rain_prob',ir,slab_size) des_rain_prob= get_sds(sd_id,'des_rain_prob',ir,slab_size) asc_rain_flag= get_sds(sd_id,'asc_rain_flag',ir,slab_size) des_rain_flag= get_sds(sd_id,'des_rain_flag',ir,slab_size) ; Select latitude and longitudes read,'Enter the minimum and maximum latitudes [-90,90]:',lat1,lat2 if ((lat1 lt -90) or (lat2 lt -90) or (lat1 gt 90) or $ (lat2 gt 90)) then begin print,'ERROR: Latitudes must be between -90 and 90' stop endif ; Make sure that lat2 is greater than lat1 if (lat1 gt lat2) then begin itmp=lat1 lat1=lat2 lat2=itmp endif ; The last grid point is in cell 719. Reduce lat 90. to 89.9 if (lat2 eq 90.) then begin lat2=89.9 endif read,'Enter the minimum and maximum longitudes [0,360]:',lon1,lon2 if ((lon1 lt 0) or (lon2 lt 0) or (lon1 gt 360) or $ (lon2 gt 360)) then begin print,'ERROR: Longitudes must be between 0 and 360' stop endif ; Make sure that lon2 is greater than lon1 if (lon1 gt lon2) then begin itmp=lon1 lon1=lon2 lon2=itmp endif ; The last grid point is in cell 1439. Wrapping is not done here, ; so 360, must be reduced to 359.9. if (lon2 eq 360.) then begin lon2=359.9 endif ; Determine grid points from the latitudes and longitudes dx=(360./1440.) ii1=fix(lon1/dx) ii2=fix(lon2/dx) dy=(180./720.) jj1=fix((lat1+90.)/dy) jj2=fix((lat2+90.)/dy) ; Print results to screen print,'ASCENDING PASS (DAYTIME)' print,' LON LAT SPD U V SPD2 COUNT TIME RAIN PROB RAIN FLAG' for i=ii1,ii2 do begin for j=jj1,jj2 do begin ; Calculate the center lat and lon of the grid box lon=i*dx+(dx/2) lat=((j*dy)-90.)+(dy/2) if (asc_wvc_count(i,j) gt 0) then begin print,lon,lat,asc_avg_wind_speed(i,j),asc_avg_wind_vel_u(i,j), $ asc_avg_wind_vel_v(i,j),asc_avg_wind_speed_sq(i,j), $ asc_wvc_count(i,j),asc_time_frac(i,j),asc_rain_prob(i,j), $ asc_rain_flag(i,j), $ format='(f9.5,x,f9.5,x,f7.2,x,f7.2,x,f7.2,x,f7.2,2x,f2.0,f9.5,x,f7.3,7x,f2.0)' endif endfor endfor print,' ' print,'DESCENDING PASS (NIGHTTIME)' print,' LON LAT SPD U V SPD2 COUNT TIME RAIN PROB RAIN FLAG' for i=ii1,ii2 do begin for j=jj1,jj2 do begin ; Calculate the center lat and lon of the grid box lon=i*dy+(dx/2) lat=((j*dy)-90.)+(dy/2) if (des_wvc_count(i,j) gt 0) then begin print,lon,lat,des_avg_wind_speed(i,j),des_avg_wind_vel_u(i,j), $ des_avg_wind_vel_v(i,j),des_avg_wind_speed_sq(i,j), $ des_wvc_count(i,j),des_time_frac(i,j),des_rain_prob(i,j), $ des_rain_flag(i,j), $ format='(f9.5,x,f9.5,x,f7.2,x,f7.2,x,f7.2,x,f7.2,2x,f2.0,f9.5,x,f7.3,7x,f2.0)' endif endfor endfor end FUNCTION ls_qscat,yr=yr,doy=doy,dr=dr lat = [41,49] lon = [-93,-76]+360 dx=(360./1440.) ii = fix(lon/dx) dy=(180./720.) jj=fix((lat+90.)/dy) IF n_elements(yr) EQ 0 THEN yr = 2006 IF n_elements(doy) EQ 0 THEN doy = 193 IF n_elements(dr) EQ 0 THEN dr = mydocs()+'superior/quikscat/raw/' fname = file_search(dr+'QS_XWGRD3_'+string(yr,format='(i4.4)')+string(doy,format='(i3.3)')+'*.gz',count=c) IF c GT 0 THEN BEGIN filename = fname[0] newfilename = strmid(filename,0,strlen(filename)-3) spawn,'gunzip -c '+filename+' > '+newfilename sd_id=HDF_SD_START(newfilename,/READ) index=HDF_SD_NAMETOINDEX(sd_id,'asc_wvc_count') sds_id=HDF_SD_SELECT(sd_id,index) HDF_SD_GETINFO,sds_id,ndims=rank,dims=dims,label=name, $ type=data_type,caldata=cal slab_size=strtrim(dims(1),2) ir=0 asc_avg_wind_speed= get_sds(sd_id,'asc_avg_wind_speed',ir,slab_size) des_avg_wind_speed= get_sds(sd_id,'des_avg_wind_speed',ir,slab_size) asc_avg_wind_vel_u= get_sds(sd_id,'asc_avg_wind_vel_u',ir,slab_size) des_avg_wind_vel_u= get_sds(sd_id,'des_avg_wind_vel_u',ir,slab_size) asc_avg_wind_vel_v= get_sds(sd_id,'asc_avg_wind_vel_v',ir,slab_size) des_avg_wind_vel_v= get_sds(sd_id,'des_avg_wind_vel_v',ir,slab_size) asc_avg_wind_speed_sq= get_sds(sd_id,'asc_avg_wind_speed_sq',ir,slab_size) des_avg_wind_speed_sq= get_sds(sd_id,'des_avg_wind_speed_sq',ir,slab_size) asc_wvc_count= get_sds(sd_id,'asc_wvc_count',ir,slab_size) des_wvc_count= get_sds(sd_id,'des_wvc_count',ir,slab_size) asc_time_frac= get_sds(sd_id,'asc_time_frac',ir,slab_size) des_time_frac= get_sds(sd_id,'des_time_frac',ir,slab_size) asc_rain_prob= get_sds(sd_id,'asc_rain_prob',ir,slab_size) des_rain_prob= get_sds(sd_id,'des_rain_prob',ir,slab_size) asc_rain_flag= get_sds(sd_id,'asc_rain_flag',ir,slab_size) des_rain_flag= get_sds(sd_id,'des_rain_flag',ir,slab_size) hdf_sd_end, sd_id file_delete,newfilename c1 = asc_wvc_count[ii[0]:ii[1],jj[0]:jj[1]] c2 = des_wvc_count[ii[0]:ii[1],jj[0]:jj[1]] u1 = asc_avg_wind_vel_u[ii[0]:ii[1],jj[0]:jj[1]] u2 = des_avg_wind_vel_u[ii[0]:ii[1],jj[0]:jj[1]] v1 = asc_avg_wind_vel_v[ii[0]:ii[1],jj[0]:jj[1]] v2 = des_avg_wind_vel_v[ii[0]:ii[1],jj[0]:jj[1]] s1 = asc_avg_wind_speed[ii[0]:ii[1],jj[0]:jj[1]] s2 = des_avg_wind_speed[ii[0]:ii[1],jj[0]:jj[1]] sq1 = asc_avg_wind_speed_sq[ii[0]:ii[1],jj[0]:jj[1]] sq2 = des_avg_wind_speed_sq[ii[0]:ii[1],jj[0]:jj[1]] tm1 = asc_time_frac[ii[0]:ii[1],jj[0]:jj[1]] tm2 = des_time_frac[ii[0]:ii[1],jj[0]:jj[1]] r1 = asc_rain_prob[ii[0]:ii[1],jj[0]:jj[1]] r2 = des_rain_prob[ii[0]:ii[1],jj[0]:jj[1]] f1 = asc_rain_flag[ii[0]:ii[1],jj[0]:jj[1]] f2 = des_rain_flag[ii[0]:ii[1],jj[0]:jj[1]] bval = where(c1 EQ 0,nbval) IF nbval GT 0 THEN BEGIN u1[bval] = nan() v1[bval] = nan() s1[bval] = nan() sq1[bval] = nan() tm1[bval] = nan() r1[bval] = nan() f1[bval] = nan() ENDIF bval2 = where(c2 EQ 0,nbval2) IF nbval2 GT 0 THEN BEGIN u2[bval2] = nan() v2[bval2] = nan() s2[bval2] = nan() sq2[bval2] = nan() tm2[bval2] = nan() r2[bval2] = nan() f2[bval2] = nan() ENDIF lons = numgen(ii[0]*dx+(dx/2),ii[1]*dx+(dx/2),dx)-360. lats = numgen( ((jj[0]*dy)-90.)+(dy/2), ((jj[1]*dy)-90.)+(dy/2),dy) lonarr = u1 FOR i = 0,n_elements(u1[0,*])-1 DO lonarr[*,i] = lons latarr = u1 FOR i = 0,n_elements(u1[*,0])-1 DO latarr[i,*] = lats tmval = where(finite(tm2) AND (tm2 LT 0.25),ntmv) IF ntmv GT 0 THEN tm2[tmval]++ tm1m = mean(tm1,/nan) tm2m = mean(tm2,/nan) output = {good:1l,u1:u1,u2:u2,v1:v1,v2:v2,s1:s1,s2:s2,sq1:sq1,sq2:Sq2,tm1:tm1,tm2:tm2,tm1m:tm1m,tm2m:tm2m,r1:r1,r2:r2,f1:f1,f2:f2,lonarr:lonarr,latarr:latarr} ENDIF ELSE BEGIN print,' File '+fname+' not found' output = {good:0l} ENDELSE ;stop return,output END PRO conv_ls_Qscat ;arrays ;netcdf? first sav file ;yr = 1999,2008 ;doy = 1999-2008 * 365 ex 366 in 00 and 04 * 2 ;tm = as many elems as doy ;lonarr,latarr = 69,33 ;u,v,s,sq,r,f = 69,33,elems doy = findgen(366)+1 yr = findgen(10)+1999 tm = make_array(10,366,2,/float,value=nan()) ; yrtm = tm u = make_array(69,33,10,366,2,/float,value=nan()) v = u s = u sq = u r = u f = u time = u ; print,'Reading file' ; restore,mydocs()+'superior/quikscat/greatlake_qs.sav' ; yrtm = tm FOR y = 0,9 DO BEGIN ; FOR y = 9,9 DO BEGIN IF (y-1) MOD 4 EQ 0 THEN diy = 366.d ELSE diy = 365.d FOR d = 0,365 DO BEGIN ; FOR d = 243,365 DO BEGIN print,'On Year ',y+1999,' Day ',d+1 ; output = ls_qscat(yr=y+1999,doy=d+1,dr='/Volumes/MyBook_Mac/Quikscat/') output = ls_qscat(yr=y+1999,doy=d+1) IF output.good EQ 1 THEN BEGIN u[*,*,y,d,0] = output.u1 v[*,*,y,d,0] = output.v1 s[*,*,y,d,0] = output.s1 sq[*,*,y,d,0] = output.sq1 r[*,*,y,d,0] = output.r1 f[*,*,y,d,0] = output.f1 u[*,*,y,d,1] = output.u2 v[*,*,y,d,1] = output.v2 s[*,*,y,d,1] = output.s2 sq[*,*,y,d,1] = output.sq2 r[*,*,y,d,1] = output.r2 f[*,*,y,d,1] = output.f2 time[*,*,y,d,0] = output.tm1 time[*,*,y,d,1] = output.tm2 tm[y,d,0] = output.tm1m*24. tm[y,d,1] = output.tm2m*24. ; yrtm[y,d,0] = (y+1999.d)+((d+1)/diy)+(output.tm1m/diy) ; yrtm[y,d,1] = (y+1999.d)+((d+1)/diy)+(output.tm2m/diy) IF n_elements(lonarr) EQ 0 THEN lonarr = output.lonarr IF n_elements(latarr) EQ 0 THEN latarr = output.latarr ENDIF ENDFOR ENDFOR ;stop save,filename=mydocs()+'superior/quikscat/greatlake_qs2.sav',u,v,s,sq,r,f,tm,latarr,lonarr,doy,yr,time ;for each doy/tm call ls_qscat ;if output.good eq 1 then copy vals ;save save file ;get mean for lake sup only for each point, plot it up ;write to netcdf eventually END PRO ls_qscatfixtime restore,mydocs()+'superior/quikscat/greatlake_qs.sav' diff = [0.0173841,0.0174711,0.0176535,0.01789791,0.0179126,0.0179047,0.0179359,0.0179189,0.0179687,0.0179639] FOR y = 0,9 DO BEGIN ttime = time[*,*,y,*,1] btime = where(ttime GE 1.0,nbt) IF nbt GT 0 THEN ttime[btime] -= diff[y] time[*,*,y,*,1] = ttime ENDFOR ;fix mean time1 ntime = reform(time[*,*,*,*,1],69.*33.,10.,366.) mntime_n = total(ntime,1,/nan)/total(finite(ntime),1,/nan) dtime = reform(time[*,*,*,*,0],69.*33.,10.,366.) mntime_d = total(dtime,1,/nan)/total(finite(dtime),1,/nan) tm[*,*,0] = mntime_d*24. tm[*,*,1] = mntime_n*24. save,filename=mydocs()+'superior/quikscat/greatlake_qs2.sav',u,v,s,sq,r,f,tm,latarr,lonarr,doy,yr,time stop END PRO ls_qscattime,time=time IF n_elements(time) EQ 0 THEN restore,mydocs()+'superior/quikscat/greatlake_qs.sav' ;diffday = [0.0] ;diff = [0.0] FOR y = 0,9 DO BEGIN FOR d = 0,365 DO BEGIN ntime = reform(time[*,*,y,d,1]) btime = where(ntime GE 1.0,nbt) gtime = where(ntime GT 0.5 AND ntime LT 1.0,ngt) IF (nbt GT 10) AND (ngt GT 10) THEN BEGIN IF n_elements(ntime1) EQ 0 THEN ntime1 = mean(ntime,/nan) ELSE ntime1 = [ntime1,mean(ntime,/nan)] IF n_elements(ntime2) EQ 0 THEN ntime2 = mean(ntime[btime],/nan) ELSE ntime2 = [ntime2,mean(ntime[btime],/nan)] IF n_elements(ntime3) EQ 0 THEN ntime3 = mean(ntime[gtime]) ELSE ntime3 = [ntime3,mean(ntime[gtime],/nan)] ktime = ntime ktime[btime] -= 0.0179189 IF n_elements(ntime4) EQ 0 THEN ntime4 = mean(ktime,/nan) ELSE ntime4 = [ntime4,mean(ktime,/nan)] tdiff = mean(ntime[btime],/nan)-mean(ntime[gtime],/nan) IF n_elements(diffday) EQ 0 THEN diffday = d ELSE diffday = [diffday,d] IF n_elements(diffyr) EQ 0 THEN diffyr = y ELSE diffyr = [diffyr,y] IF n_elements(diff) EQ 0 THEN diff = tdiff ELSE diff = [diff,tdiff] ENDIF ENDFOR ENDFOR stop END PRO ls_qscatcdf,yrs=yrs ;restore the file restore,mydocs()+'superior/quikscat/greatlake_qs.sav' IF n_elements(yrs) EQ 0 THEN yrs = [2008] FOR theyr = 0,n_elements(yrs)-1 DO BEGIN yr = yrs[theyr] print,'On Year ',yr yrloc = yr-1999l yrstr = string(yr,format='(i4.4)') basetimes = [915148800l, 946684800l ,978307200l, 1009843200l,1041379200l ,1072915200l ,1104537600l,1136073600l,1167609600l,1199145600l] dd = congrid(doy,366.*2,/center)-1 newtm = ((dd + (reform(transpose(reform(tm[yrloc,*,*])),732)/24.))*86400.0) newtm2 = newtm badtm = where((newtm[1:*:2] LT newtm[0:*:2]) AND finite(newtm[1:*:2]) AND finite(newtm[0:*:2]),nbt) IF nbt GT 0 THEN newtm[(badtm*2)+1]+=86400. gtm = where(finite(newtm),ntm) newtm = long(newtm[gtm]) newu = (reform(transpose(reform(u[*,*,yrloc,*,*]),[0,1,3,2]),69,33,732))[*,*,gtm] newv = (reform(transpose(reform(v[*,*,yrloc,*,*]),[0,1,3,2]),69,33,732))[*,*,gtm] newws = (reform(transpose(reform(s[*,*,yrloc,*,*]),[0,1,3,2]),69,33,732))[*,*,gtm] newsq = (reform(transpose(reform(sq[*,*,yrloc,*,*]),[0,1,3,2]),69,33,732))[*,*,gtm] newr = (reform(transpose(reform(r[*,*,yrloc,*,*]),[0,1,3,2]),69,33,732))[*,*,gtm] newf = (reform(transpose(reform(f[*,*,yrloc,*,*]),[0,1,3,2]),69,33,732))[*,*,gtm] newtime = (reform(transpose(reform(time[*,*,yrloc,*,*]),[0,1,3,2]),69,33,732))[*,*,gtm] FOR i = 0,32 DO BEGIN FOR j = 0,68 DO BEGIN newtime[j,i,*]+=(dd+1) ENDFOR ENDFOR bval = where(~finite(newu),nbval) & IF nbval GT 0 THEN newu[bval]=-999 bval = where(~finite(newv),nbval) & IF nbval GT 0 THEN newv[bval]=-999 bval = where(~finite(newws),nbval) & IF nbval GT 0 THEN newws[bval]=-999 bval = where(~finite(newsq),nbval) & IF nbval GT 0 THEN newsq[bval]=-999 bval = where(~finite(newr),nbval) & IF nbval GT 0 THEN newr[bval]=-999 bval = where(~finite(newf),nbval) & IF nbval GT 0 THEN newf[bval]=-999 bval = where(~finite(newtime),nbval) & IF nbval GT 0 THEN newf[bval]=-999 ; id = ncdf_open('/data/projects/superior/quikscat/glake_qscat_2008.nc') ; ncdf_varget,id,'time_offset',toff id2 = ncdf_create('/data/projects/superior/quikscat/glake_qscat_'+yrstr+'.nc',/clobber) ncdf_attput,id2,'TITLE','Great Lakes Quikscat Level 3 Wind Data',/global ncdf_attput,id2,'YEAR',yr,/global ncdf_attput,id2,'BASE_TIME_DATE',yrstr+'-01-01 00:00:00 UTC',/global ncdf_attput,id2,'START_DATE',yrstr+'-01-01 00:00:00 UTC',/global ncdf_attput,id2,'END_DATE',yrstr+'-01-01 00:00:00 UTC',/global ncdf_attput,id2,'LON_GRID_DIMENSION',69,/global ncdf_attput,id2,'LAT_GRID_DIMENSION',33,/global ncdf_attput,id2,'TIME_DIMENSION',ntm,/global ncdf_attput,id2,'MISSING_DATA',-999,/global ncdf_attput,id2,'CREATOR','Ankur Desai, UW-Madison, desai@aos.wisc.edu',/global lonid = ncdf_dimdef(id2,'lon',69) latid = ncdf_dimdef(id2,'lat',33) tmid = ncdf_dimdef(id2,'time',ntm) btid = ncdf_vardef(id2,'basetime',/long) time_offset = ncdf_vardef(id2,'time_offset',tmid,/long) long = ncdf_vardef(id2,'longitude',lonid,/float) lati = ncdf_vardef(id2,'latitude',latid,/float) u_var = ncdf_vardef(id2,'U',[tmid,lonid,latid],/float) v_var = ncdf_vardef(id2,'V',[tmid,lonid,latid],/float) tm_var = ncdf_vardef(id2,'DAY_FRAC',[tmid,lonid,latid],/float) ws_var = ncdf_vardef(id2,'WIND_SPEED',[tmid,lonid,latid],/float) wss_var = ncdf_vardef(id2,'WIND_SPEED_SQ',[tmid,lonid,latid],/float) rp_var = ncdf_vardef(id2,'RAIN_PROB',[tmid,lonid,latid],/float) rf_var = ncdf_vardef(id2,'RAIN_FLAG',[tmid,lonid,latid],/float) ncdf_control,id2,/endef ncdf_varput,id2,btid,basetimes[yrloc] ncdf_varput,id2,time_offset,reform(newtm) ncdf_varput,id2,long,reform(lonarr[*,0]) ncdf_varput,id2,lati,reform(latarr[0,*]) ncdf_varput,id2,u_var,transpose(reform(newu),[2,0,1]) ncdf_varput,id2,v_var,transpose(reform(newv),[2,0,1]) ncdf_varput,id2,tm_var,transpose(reform(newtime),[2,0,1]) ncdf_varput,id2,ws_var,transpose(reform(newws),[2,0,1]) ncdf_varput,id2,wss_var,transpose(reform(newsq),[2,0,1]) ncdf_varput,id2,rp_var,transpose(reform(newr),[2,0,1]) ncdf_varput,id2,rf_var,transpose(reform(newf),[2,0,1]) ncdf_close,id2 ENDFOR END ;id=ndcf_create('greatlake_qs.nc') ;ncdf_attput,id,'Title','Hello',/global ;... missing_data = -999 ;ncdf_control,id,/fill ;xid = ncdf_dimdef(id,'x',100) ;... ;vid = ncdf_vardef(id,'varname',[yid,xid],/byte) /float ... ;ncdf_attput,id,vid,'units','m/s' ;ncdf_control,id,/endef ;ncdf_varput,id,vid,var ;ncdf_close,id ;base_time in seconds since 1970-1-1 00 (variable) ;time in seconds since base_time ;latitude (y) ;longitude (x) PRO ls_testwfunc,deltaT=deltaT,h=h COMMON sharew, dT, z IF n_elements(deltaT) NE 0 THEN BEGIN dT = deltaT IF n_elements(h) EQ 0 THEN z = 10 ELSE z = h Urun = ls_wfunc(1.0) Uguess = 1.0 print,newton(Uguess,'ls_wfunc',/double) ENDIF ELSE BEGIN z = 3 ; dT = -10 ; print,ls_wfunc(30) Uarr = fltarr(20,10) dTs = (findgen(20)-20)/8.0 FOR zs = 1,10 DO BEGIN print,zs z = zs FOR x=0,19 DO BEGIN dT = dTs[x] Urun = ls_wfunc(1.0) Uguess = 1.0 Uarr[x,zs-1] = newton(Uguess,'ls_wfunc',/double) ; stop ENDFOR ENDFOR stop ENDELSE END FUNCTION ls_wfunc,x COMMON sharew U = (x[0]>0.001) ; print,dT,z ; IF n_elements(dT) EQ 0 THEN dT = -5 ; IF n_elements(z) EQ 0 THEN z = 10 Cs = 4.4e-4 * U^0.55 ; Cs = 0.86e-3 ; ustar = 0.0409 * U + 0.0946 ustar = sqrt(4.4e-4 * (U^2.55)) z0 = (0.015/9.8) * ustar^2 ; z0 = 0.016 ; H = Cs * U * dT H = 2.6e-3 + 0.864e-3 * U * dT kgT = 0.013 ; L = (-1.0) * ustar^3 / (kgT * H) ; psi = (4.7) * z / L psi = (-4.7) * z * kgT * H / (ustar^3) outval = ( (ustar/0.4) * (alog(z/z0) + psi) ) - U return,outval END FUNCTION ls_windprofile ;M-O theory ; ustar = 0.1 kgT = 0.013 ;U = (u*/k)(ln(z/z0)-psi(z/L)) ;L = -u*^3 / ( (kg/To * (H0/rhocp))) TlTa = findgen(10)-5 ;psi(z/L) = -5z/l ;U = 4.4e-4/0.4 * U^2.55 * [ln(z/z0)+5z/L] ;Qs = Cs U deltaT ;z0 = 0.015 u*^2 / g ;u* = 4.4e-4 * U^2.55 END