; Program to read the 2011 HS3 HAMSR and AVAPS data ; and compare the two dataets ; Andrew Kren ; March 11, 2016 begin ; cases case = "20110909" path_avaps = "/scratch4/BMC/shout/Andrew.Kren/avaps_data/2011/Pacific/" ;+ case + "/" ; dropsonde path path_hamsr = "/scratch4/BMC/shout/Andrew.Kren/hamsr_data/2011/HS3/L2/" ;+ case + "/" ; hamsr path ; list files in directory to determine how many dropsondes to read fils = systemfunc("csh -c 'cd " + path_avaps + " ; ls *.nc'") fils2 = systemfunc("ls " + path_hamsr + "*.nc") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; end of user's settings ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function closest_val_AnyOrder(xVal[1]:numeric, x:numeric) local xAbsDif, xMinVal, iClose begin xAbsDif = abs(xVal-x) iClose = minind(xAbsDif) return(iClose) ; "first occurrence" end ; extract dates from both datasets year_avaps = tointeger(str_get_cols(fils,1,4)) ; year, month, day, hour, minute, sec in UTC, same for hamsr month_avaps = tointeger(str_get_cols(fils,5,6)) day_avaps = tointeger(str_get_cols(fils,7,8)) hour_avaps = tointeger(str_get_cols(fils,10,11)) minute_avaps = tointeger(str_get_cols(fils,12,13)) sec_avaps = tointeger(str_get_cols(fils,14,15)) year_avaps_file = str_get_cols(fils,1,4) ; year, month, day, hour, minute, sec in UTC, same for hamsr month_avaps_file = str_get_cols(fils,5,6) day_avaps_file = str_get_cols(fils,7,8) hour_avaps_file = str_get_cols(fils,10,11) minute_avaps_file = str_get_cols(fils,12,13) time_avaps_plot = str_get_cols(fils,10,13) print(fils2) ; read in all HAMSR data q = addfiles(fils2,"r") factor = 0.001 lat_hamsr = factor*(q[:]->lat) ; all lats, in (time,lat) format lon_hamsr = factor*(q[:]->lon) factor = 0.1 lev_hamsr = factor*(q[:]->ham_pres_levels) ; mb time_hamsr = q[:]->time ; all times time_hamsr@units = "seconds since 2000-01-01 00:00:00.0" temp_hamsr = factor*(short2flt(q[:]->ham_airT))-273.15 ; degC exit rh_hamsr = factor*(q[:]->ham_airRH) ; % factor = 0.001 abshum_hamsr = factor*(q[:]->ham_airQ) ; g/m3 factor = 0.01 incidence = factor*(q[:]->inc) ; global hawk incidence angle exit do i=0,dimsizes(fils)-1 ; loop over all dropsondes for this case and compare to hamsr print("Reading AVAPS file: " + fils(i)) ; find hamsr file that closely matches the avaps file closest_index = closest_val_AnyOrder( hour_avaps(i), hour_hamsr_1) print("Closest HAMSR file to AVAPS: " + fils2(closest_index)) ; read data ; read file and save to array of (lat,lon) or (lev,lat,lon) w = addfile(path_avaps+fils(i),"r") lat_avaps = w->lat lon_avaps = w->lon lev_avaps = w->pres ; mb time_avaps = w->time hr_avaps = w->hr min_avaps = w->minute second_avaps = w->sec temp_avaps = w->temp ; degC rh_avaps = w->rh ; % dew_avaps = w->dew ; degC ; calculate absolute humidity to compare with HAMSR absolute humidity vapor_pres = 6.11*10.0^((7.5*dew_avaps)/(237.3 + dew_avaps))*100.0 ; Pascals temp_kelvin = temp_avaps + 273.15 ; temperature in degK r_w = 461.5 ; J/kg*K abs_hum = vapor_pres / (temp_kelvin * r_w) * 1000.0 ; g/m3 q = addfile(path_hamsr+fils2(closest_index),"r") lat_hamsr = q->lat lat_2 = tofloat(lat_hamsr) lat_hamsr_1 = lat_2 * 0.001 ; pixel lat/lon lon_hamsr = q->lon lon_2 = tofloat(lon_hamsr) lon_hamsr_1 = lon_2 * 0.001 lev_hamsr = q->ham_pres_levels ; mb lev_2 = tofloat(lev_hamsr) lev_hamsr_1 = lev_2 * 0.1 time_hamsr = q->time time_hamsr@units = "seconds since 2000-01-01 00:00:00.0" temp_hamsr = q->ham_airT temp_2 = tofloat(temp_hamsr) temp_hamsr_1 = (temp_2 * 0.1) - 273.15 ; degC rh_hamsr = q->ham_airRH ; % rh_2 = tofloat(rh_hamsr) rh_hamsr_1 = rh_2 * 0.01 abshum_ham = q->ham_airQ abshum_ham2 = tofloat(abshum_ham) abshum_ham1 = abshum_ham2 * 0.001 ; g/m3 incidence = q->inc ; Global Hawk incidence angle incidence_2 = tofloat(incidence) incidence_1 = incidence_2 * 0.01 ; should be 0.01 factor per discussion with Gary Wick ; print(incidence_1(0,:)) time_hamsr_2 = cd_calendar(time_hamsr,0) hour_hamsr = time_hamsr_2(:,3) minute_hamsr = time_hamsr_2(:,4) ; find hamsr minute cloest to time of release of dropsonde from aircraft, plotting that closest minute to the end closest_index_min = closest_val_AnyOrder( min_avaps(0), minute_hamsr) print("Closest HAMSR minute to AVAPS: " + minute_hamsr(closest_index_min)) ; filter out -999 values and other things... for plotting indices_lev = ind(lev_avaps.ne.-999) lev_avaps_2 = lev_avaps(indices_lev) temp_avaps_2 = temp_avaps(indices_lev) rh_avaps_2 = rh_avaps(indices_lev) rh_avaps_2@long_name = "Relative Humidity (%)" lev_avaps_2@long_name = "Pressure (hPa)" temp_avaps_2@long_name = "Temperature (degC)" abs_hum@long_name = "Absolute Humidity (g/m3)" abs_hum_2 = abs_hum(indices_lev) temp_hamsr_1 = where(temp_hamsr_1.lt.-300,-999,temp_hamsr_1) temp_hamsr_1@_FillValue = -999 rh_hamsr_1 = where(rh_hamsr_1.lt.0,-999,rh_hamsr_1) rh_hamsr_1@_FillValue = -999 abshum_ham1 = where(abshum_ham1.lt.0,-999,abshum_ham1) abshum_ham1@_FillValue = -999 temp_hamsr_1 = temp_hamsr_1(:,:,::-1) rh_hamsr_1 = rh_hamsr_1(:,:,::-1) ; reverse levels since goes from 1000 to 0.1 lev_hamsr_1 = lev_hamsr_1(::-1) abshum_ham1 = abshum_ham1(:,:,::-1) temp_hamsr_2 = dim_avg_n_Wrap(temp_hamsr_1,1) ; time x lev array, average over lat rh_hamsr_2 = dim_avg_n_Wrap(rh_hamsr_1,1) ; time x lev array, average over lat printVarSummary(rh_hamsr_2) print(lev_hamsr_1) exit ; find acceptable incidence angles where angle +/- 3 degrees to get (near) nadir retrievals from HAMSR do k=0,dimsizes(time_hamsr)-1 acceptable_angles = ind(incidence_1(k,:).gt.3.0) acceptable_angles@_FillValue = -999 temp_hamsr_1(k,acceptable_angles,:) = -999 rh_hamsr_1(k,acceptable_angles,:) = -999 abshum_ham1(k,acceptable_angles,:) = -999 delete(acceptable_angles) end do temp_hamsr_2 = dim_avg_n_Wrap(temp_hamsr_1,1) rh_hamsr_2 = dim_avg_n_Wrap(rh_hamsr_1,1) abshum_ham_2 = dim_avg_n_Wrap(abshum_ham1,1) ; absolute humidity ; interpolate the avaps data to the pressure level of HAMSR linlog = 2 ; ln(p) interpolation temp_avaps_2_int = int2p(lev_avaps_2,temp_avaps_2,lev_hamsr_1,linlog) abs_hum_2_int = int2p(lev_avaps_2,abs_hum_2,lev_hamsr_1,linlog) rh_avaps_2_int = int2p(lev_avaps_2,rh_avaps_2,lev_hamsr_1,linlog) ; plot vertical profiles of temperature and relative humidity, with HAMSR and AVAPS overlayed ; plot temp and RH on separate profiles ;************************************************ ; create plots, Temperature, then RH ;************************************************ wks = gsn_open_wks ("png","avaps_hamsr_temp_"+year_avaps_file(i)+month_avaps_file(i)+day_avaps_file(i)+hour_avaps_file(i)+minute_avaps_file(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "AVAPS and HAMSR " + month_avaps(i) + "/" + day_avaps(i) + "/" + year_avaps(i) + " " + time_avaps_plot(i) + " UTC" res@trYReverse = True ; reverse Y-axis res@gsnFrame = False ; don't advance frame yet colors_plot = (/"red","blue","green"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = -80 res@trXMaxF = 30 ; res@trYLog = True res@trYMinF = 0.1 res@trYMaxF = 1000.0 res@tmYLMode = "Explicit" ; Force labels where we want them res@tmYLValues = (/10,100,200,300,400,500,600,700,800,900,1000/) res@tmYLLabels = "" + res@tmYLValues res@tiYAxisString = "Pressure (hPa)" res@tiXAxisString = "Temperature (degC)" res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,temp_avaps_2_int,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 plot1 = gsn_csm_xy(wks,temp_hamsr_2(closest_index_min,:),lev_hamsr_1,res) ; Attach a legend lgres = True lgres@lgLineColors = colors_plot lgres@lgItemType = "Lines" ; show lines only (default) lgres@lgLabelFontHeightF = .08 ; legend label font thickness lgres@vpWidthF = 0.13 ; width of legend (NDC) lgres@vpHeightF = 0.10 ; height of legend (NDC) lgres@lgPerimThicknessF = 2.2 ; thicken the box perimeter lgres@lgMonoDashIndex = True lgres@lgDashIndex = 0 labels = (/"AVAPS","HAMSR"/);,"HAMSR 2"/) legend = gsn_create_legend (wks, 2, labels,lgres) ; ; Use gsn_add_annotation to attach this legend to our existing plot. ; This way, if we resize the plot, the legend will stay with the ; plot and be resized automatically. ; ; Point (0,0) is the dead center of the plot. Point (0,.5) is center, ; flush bottom. Point (0.5,0.5) is flush bottom, flush right. ; amres = True amres@amJust = "BottomRight" ; Use bottom right corner of box ; for determining its location. amres@amParallelPosF = 0.5 ; Move legend to right amres@amOrthogonalPosF = -0.3 ; Move legend down. annoid = gsn_add_annotation(plot0,legend,amres) ; add legend to plot overlay(plot0,plot1) ; plot additional hamsr retrievals after the closest minute ;do j=closest_index_min+1,dimsizes(minute_hamsr)-1,10 ; res@xyLineColor = colors_plot(2) ; plot2 = gsn_csm_xy(wks,temp_hamsr_2(j,:),lev_hamsr_1,res) ; overlay(plot0,plot2) ;end do res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot3 = gsn_csm_xy(wks,temp_avaps_2_int,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) plot4 = gsn_csm_xy(wks,temp_hamsr_2(closest_index_min,:),lev_hamsr_1,res) overlay(plot0,plot3) overlay(plot0,plot4) draw(plot0) frame(wks) ;;;;;;;;;;;;;;;;;;; ; relative humidity ;;;;;;;;;;;;;;;;;;; wks = gsn_open_wks ("png","avaps_hamsr_rh_"+year_avaps_file(i)+month_avaps_file(i)+day_avaps_file(i)+hour_avaps_file(i)+minute_avaps_file(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "AVAPS and HAMSR " + month_avaps(i) + "/" + day_avaps(i) + "/" + year_avaps(i) + " " + time_avaps_plot(i) + " UTC" res@trYReverse = True ; reverse Y-axis res@gsnFrame = False ; don't advance frame yet colors_plot = (/"red","blue","green"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = 100 ; res@trYLog = True res@trYMinF = 0.1 res@trYMaxF = 1000.0 res@tmYLMode = "Explicit" ; Force labels where we want them res@tmYLValues = (/10,100,200,300,400,500,600,700,800,900,1000/) res@tmYLLabels = "" + res@tmYLValues res@tiYAxisString = "Pressure (hPa)" res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) res@tiXAxisString = "Relative Humidity (%)" plot0 = gsn_csm_xy(wks,rh_avaps_2_int,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 plot1 = gsn_csm_xy(wks,rh_hamsr_2(closest_index_min,:),lev_hamsr_1,res) ; Attach a legend lgres = True lgres@lgLineColors = colors_plot lgres@lgItemType = "Lines" ; show lines only (default) lgres@lgLabelFontHeightF = .08 ; legend label font thickness lgres@vpWidthF = 0.13 ; width of legend (NDC) lgres@vpHeightF = 0.10 ; height of legend (NDC) lgres@lgPerimThicknessF = 2.2 ; thicken the box perimeter lgres@lgMonoDashIndex = True lgres@lgDashIndex = 0 labels = (/"AVAPS","HAMSR"/);,"HAMSR 2"/) legend = gsn_create_legend (wks, 2, labels,lgres) ; ; Use gsn_add_annotation to attach this legend to our existing plot. ; This way, if we resize the plot, the legend will stay with the ; plot and be resized automatically. ; ; Point (0,0) is the dead center of the plot. Point (0,.5) is center, ; flush bottom. Point (0.5,0.5) is flush bottom, flush right. ; amres = True amres@amJust = "BottomRight" ; Use bottom right corner of box ; for determining its location. amres@amParallelPosF = 0.5 ; Move legend to right amres@amOrthogonalPosF = -0.3 ; Move legend down. annoid = gsn_add_annotation(plot0,legend,amres) ; add legend to plot overlay(plot0,plot1) ; plot additional hamsr retrievals after the closest minute ;do j=closest_index_min+1,dimsizes(minute_hamsr)-1,10 ; res@xyLineColor = colors_plot(2) ; plot2 = gsn_csm_xy(wks,rh_hamsr_2(j,:),lev_hamsr_1,res) ; overlay(plot0,plot2) ;end do res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot3 = gsn_csm_xy(wks,rh_avaps_2_int,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) plot4 = gsn_csm_xy(wks,rh_hamsr_2(closest_index_min,:),lev_hamsr_1,res) overlay(plot0,plot3) overlay(plot0,plot4) draw(plot0) frame(wks) delete(res) delete(plot0) delete(plot1) ;delete(plot2) delete(plot3) delete(plot4) ;;;;;;;;;;;; ; absolute humidity ;;;;;;;;;;;;; wks = gsn_open_wks ("png","avaps_hamsr_abshum_"+year_avaps_file(i)+month_avaps_file(i)+day_avaps_file(i)+hour_avaps_file(i)+minute_avaps_file(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "AVAPS and HAMSR " + month_avaps(i) + "/" + day_avaps(i) + "/" + year_avaps(i) + " " + time_avaps_plot(i) + " UTC" res@trYReverse = True ; reverse Y-axis res@gsnFrame = False ; don't advance frame yet colors_plot = (/"red","blue","green"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = 15.0 ; res@trYLog = True res@trYMinF = 0.1 res@trYMaxF = 1000.0 res@tmYLMode = "Explicit" ; Force labels where we want them res@tmYLValues = (/10,100,200,300,400,500,600,700,800,900,1000/) res@tmYLLabels = "" + res@tmYLValues res@tiYAxisString = "Pressure (hPa)" res@tiXAxisString = "Absolute Humidity (g/m3)" res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,abs_hum_2_int,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 plot1 = gsn_csm_xy(wks,abshum_ham_2(closest_index_min,:),lev_hamsr_1,res) ; Attach a legend lgres = True lgres@lgLineColors = colors_plot lgres@lgItemType = "Lines" ; show lines only (default) lgres@lgLabelFontHeightF = .08 ; legend label font thickness lgres@vpWidthF = 0.13 ; width of legend (NDC) lgres@vpHeightF = 0.10 ; height of legend (NDC) lgres@lgPerimThicknessF = 2.2 ; thicken the box perimeter lgres@lgMonoDashIndex = True lgres@lgDashIndex = 0 labels = (/"AVAPS","HAMSR"/);,"HAMSR 2"/) legend = gsn_create_legend (wks, 2, labels,lgres) ; ; Use gsn_add_annotation to attach this legend to our existing plot. ; This way, if we resize the plot, the legend will stay with the ; plot and be resized automatically. ; ; Point (0,0) is the dead center of the plot. Point (0,.5) is center, ; flush bottom. Point (0.5,0.5) is flush bottom, flush right. ; amres = True amres@amJust = "BottomRight" ; Use bottom right corner of box ; for determining its location. amres@amParallelPosF = 0.5 ; Move legend to right amres@amOrthogonalPosF = -0.3 ; Move legend down. annoid = gsn_add_annotation(plot0,legend,amres) ; add legend to plot overlay(plot0,plot1) ; plot additional hamsr retrievals after the closest minute ;do j=closest_index_min+1,dimsizes(minute_hamsr)-1,10 ; res@xyLineColor = colors_plot(2) ; plot2 = gsn_csm_xy(wks,abshum_ham_2(j,:),lev_hamsr_1,res) ; overlay(plot0,plot2) ;end do res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot3 = gsn_csm_xy(wks,abs_hum_2_int,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) plot4 = gsn_csm_xy(wks,abshum_ham_2(closest_index_min,:),lev_hamsr_1,res) overlay(plot0,plot3) overlay(plot0,plot4) draw(plot0) frame(wks) ; end of plot section ; delete variables delete(w) delete(lat_avaps) delete(lon_avaps) ;delete(temp_hamsr_avg) ;delete(rh_hamsr_avg) delete(lev_avaps) delete(time_avaps) delete(hr_avaps) delete(min_avaps) delete(second_avaps) delete(temp_avaps) delete(rh_avaps) delete(dew_avaps) delete(temp_avaps_2) delete(lev_avaps_2) delete(indices_lev) delete(vapor_pres) delete(abs_hum) delete(abshum_ham) delete(abshum_ham1) delete(abshum_ham_2) delete(abshum_ham2) delete(temp_kelvin) delete(q) delete(lat_hamsr) delete(lat_2) delete(lat_hamsr_1) delete(lon_hamsr) delete(lon_2) delete(lon_hamsr_1) delete(lev_hamsr) delete(lev_2) delete(temp_hamsr_2) delete(rh_avaps_2) delete(rh_hamsr_2) delete(lev_hamsr_1) delete(time_hamsr) delete(temp_hamsr) delete(temp_2) delete(temp_hamsr_1) delete(rh_hamsr) delete(incidence) delete(incidence_1) delete(incidence_2) delete(rh_2) delete(rh_hamsr_1) delete(time_hamsr_2) delete(hour_hamsr) delete(minute_hamsr) delete(abs_hum_2) delete(abs_hum_2_int) delete(temp_avaps_2_int) end do end