; Program to read the bias-corrected HAMSR data from Shannon Brown and compare it to the dropsondes ; ; ; Andrew Kren ; May 3, 2017 begin case = (/"20161005","20161007","20161009"/) ; case names for reading data of dropsondes and hamsr, needed for directory paths ; save avaps and hasmr for plotting average profiles saved_avaps_temp = new((/150,25/),"float") saved_avaps_abshum = new((/150,25/),"float") saved_hamsr_temp = new((/150,25/),"float") saved_hamsr_abshum = new((/150,25/),"float") saved_avaps_spechum = new((/150,25/),"float") saved_hamsr_spechum = new((/150,25/),"float") ; save avaps and hasmr for plotting average profiles saved_avaps_temp@_FillValue = -999 saved_avaps_abshum@_FillValue = -999 saved_hamsr_temp@_FillValue = -999 saved_hamsr_abshum@_FillValue = -999 saved_avaps_spechum@_FillValue = -999 saved_hamsr_spechum@_FillValue = -999 gg = 0 ; counter for saving variables over all 3 storm cases do c=0,dimsizes(case)-1 ; loop over all 3 2016 cases at once print("Working on storm case: " + case(c)) path_avaps = "/scratch4/BMC/shout/Andrew.Kren/avaps_data/hrr_2016/" + case(c) + "/" ; dropsonde path path_hamsr = "/scratch4/BMC/shout/Andrew.Kren/hamsr_data/2016/2016_updated/" + case(c) + "/" ; hamsr path ; list files in directory to determine how many dropsondes to read fils = systemfunc("csh -c 'cd " + path_avaps + " ; ls *.nc'") fils2 = systemfunc("csh -c 'cd " + path_hamsr + " ; ls *.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) ; HAMSR_L2DB_realtime_20161009T0700_20161009T0800.nc ; sample name for hamsr data year_hamsr = tointeger(str_get_cols(fils2,20,23)) month_hamsr = tointeger(str_get_cols(fils2,24,25)) day_hamsr = tointeger(str_get_cols(fils2,26,27)) hour_hamsr_1 = tointeger(str_get_cols(fils2,29,30)) minute_hamsr_1 = tointeger(str_get_cols(fils2,31,32)) hour_hamsr_2 = tointeger(str_get_cols(fils2,43,44)) minute_hamsr_2 = tointeger(str_get_cols(fils2,45,46)) ; read hamsr data q = addfiles(path_hamsr+fils2,"r") lat_hamsr = q[:]->lat lat_2 = tofloat(lat_hamsr) lat_hamsr_1 = lat_2 * 0.001 ; pixel lat/lon time_hamsr = q[:]->time ; sort array across time qsort(time_hamsr) time_hamsr@units = "seconds since 2000-01-01 00:00:00.0" time_hamsr_2 = cd_calendar(time_hamsr,0) hour_hamsr = time_hamsr_2(:,3) minute_hamsr = time_hamsr_2(:,4) z = addfile(path_hamsr+fils2(0),"r") lev_hamsr = z->ham_pres_levels ; mb lev_2 = tofloat(lev_hamsr) lev_hamsr_1 = lev_2 * 0.1 temp_hamsr = q[:]->ham_airT ip = dim_pqsort(time_hamsr,2) test1 = temp_hamsr test1 = temp_hamsr(ip,:) test1&time = temp_hamsr&time(ip) temp_2 = tofloat(test1) temp_hamsr_1 = (temp_2 * 0.1) - 273.15 ; degC print(temp_hamsr_1(0:10,:)) exit 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 print(temp_hamsr_1(0:10,:)) temp_hamsr_1 = dim_pqsort_n(temp_hamsr_1,2,0) print(temp_hamsr_1(0:10,:)) exit do i=0,dimsizes(fils)-1 ; loop over all dropsondes for this case and compare to hamsr print("Reading AVAPS file: " + fils(i)) ; read in all hamsr and dropsonde data for this storm case ;;; reading dropsonde 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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; 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)) ; 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 ; reading hamsr data ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; q = addfile(path_hamsr+fils2(closest_index),"r") 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 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_start = closest_val_AnyOrder( min_avaps(0), minute_hamsr) ; closest starting time to dropsonde release temp_min_avaps = min_avaps(::-1) closest_index_end = closest_val_AnyOrder ( temp_min_avaps(0),minute_hamsr) ; closest ending time to dropsonde end time print("Closest HAMSR minute to AVAPS release: " + minute_hamsr(closest_index_start)) print("Closest HAMSR minute to AVAPS end time: " + minute_hamsr(closest_index_end)) exit difference = abs(min_avaps(0) - minute_hamsr(closest_index_min)) ; see if the difference is too far apart, greater than 10 minutes if (difference.le.25) then ; greater than 25 minutes away, toss ; 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 if (lev_hamsr_1(0).eq.1000.0) then 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) end if ; 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) ; calculate specific humidity from pressure, temperature, and relative humidity spec_hum_avaps = mixhum_ptrh(lev_avaps_2,(temp_avaps_2+273.15),rh_avaps_2,-2) ; g/kg spec_hum_avaps_int = mixhum_ptrh(lev_hamsr_1,(temp_avaps_2_int+273.15),rh_avaps_2_int,-2) ; g/kg spec_hum_avaps@_FillValue = -999 spec_hum_avaps_int@_FillValue = -999 ; now for HAMSR spec_hum_hamsr = rh_hamsr_1 do v=0,dimsizes(time_hamsr)-1 spec_hum_hamsr(v,:) = mixhum_ptrh(lev_hamsr_1,(temp_hamsr_1(v,:)+273.15),rh_hamsr_1(v,:),-2) ; g/kg end do spec_hum_hamsr@_FillValue = -999 ; 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_1(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_1(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_1(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_1(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_1(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) ;;;;;;;;;;;;;;;;;;; ; specific humidity ;;;;;;;;;;;;;;;;;;; wks = gsn_open_wks ("png","avaps_hamsr_spechum_"+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.1 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@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) res@tiXAxisString = "Specific Humidity (g/kg)" plot0 = gsn_csm_xy(wks,spec_hum_avaps_int,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 plot1 = gsn_csm_xy(wks,spec_hum_hamsr(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_1(j,:),lev_hamsr_1,res) ; overlay(plot0,plot2) ;end do res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot3 = gsn_csm_xy(wks,spec_hum_avaps_int,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) plot4 = gsn_csm_xy(wks,spec_hum_hamsr(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_ham1(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_1(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_ham1(closest_index_min,:),lev_hamsr_1,res) overlay(plot0,plot3) overlay(plot0,plot4) draw(plot0) frame(wks) ; end of plot section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; write out data to ascii text file for this time period for both HAMSR and AVAPS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; write data to output file output_text_file = "avaps_hamsr_output_"+year_avaps_file(i)+month_avaps_file(i)+day_avaps_file(i)+hour_avaps_file(i)+minute_avaps_file(i)+".txt" system("rm -rf " + output_text_file) ; header = [/"Date","Lat", "Lon", "Pres (mb)", "Drop T(K)","Ham T(K)","Drop w(g/kg)","Ham w(g/kg)","AVAPS q(g/m3)","Ham q(g/m3)","Drop RH(%)","Ham RH(%)"/] ; write_table(output_text_file, "w", header, "%12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s %12s") do o=0,dimsizes(lev_hamsr_1)-1 alist = [/year_avaps_file(i)+month_avaps_file(i)+day_avaps_file(i)+hour_avaps_file(i)+minute_avaps_file(i),lat_avaps(0),lon_avaps(0),lev_hamsr_1(o),(temp_avaps_2_int(o)+273.15),(temp_hamsr_1(closest_index_min,o)+273.15),spec_hum_avaps_int(o),spec_hum_hamsr(closest_index_min,o),abs_hum_2_int(o),abshum_ham1(closest_index_min,o),rh_avaps_2_int(o),rh_hamsr_1(closest_index_min,o)/] write_table(output_text_file, "a", alist,"%12s %12.2f %12.2f %12.2f %12.2f %12.2f %12.2f %12.2f %12.2f %12.2f %12.2f %12.2f") end do ; save varibles for this time step avaps_temp_all(gg,:) = temp_avaps_2_int avaps_abshum_all(gg,:) = abs_hum_2_int hamsr_temp_all(gg,:) = temp_hamsr_1(closest_index_min,:) hamsr_abshum_all(gg,:) = abshum_ham1(closest_index_min,:) avaps_spechum_all(gg,:) = spec_hum_avaps_int hamsr_spechum_all(gg,:) = spec_hum_hamsr(closest_index_min,:) ; delete variables delete(w) delete(lat_avaps) delete(lon_avaps) 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_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(rh_avaps_2) delete(time_hamsr) delete(temp_hamsr) delete(temp_2) delete(temp_hamsr_1) delete(rh_hamsr) delete(rh_2) delete(spec_hum_hamsr) delete(spec_hum_avaps) delete(spec_hum_avaps_int) 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) gg = gg + 1 else print("HAMSR profile too far away from AVAPS release time...skipping this dropsonde") ; delete variables delete(w) delete(abs_hum) 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_hamsr_1) delete(time_hamsr) delete(temp_hamsr) delete(temp_2) delete(temp_hamsr_1) delete(rh_hamsr) delete(rh_2) delete(vapor_pres) delete(temp_kelvin) delete(rh_hamsr_1) delete(abshum_ham) delete(abshum_ham2) delete(abshum_ham1) delete(time_hamsr_2) delete(hour_hamsr) delete(minute_hamsr) delete(lat_avaps) delete(lon_avaps) delete(lev_avaps) delete(time_avaps) delete(hr_avaps) delete(min_avaps) delete(second_avaps) delete(temp_avaps) delete(rh_avaps) delete(dew_avaps) end if end do delete(fils) delete(fils2) delete(year_avaps) delete(month_avaps) delete(day_avaps) delete(hour_avaps) delete(minute_avaps) delete(sec_avaps) delete(year_avaps_file) delete(month_avaps_file) delete(day_avaps_file) delete(hour_avaps_file) delete(minute_avaps_file) delete(time_avaps_plot) delete(year_hamsr) delete(month_hamsr) delete(day_hamsr) delete(hour_hamsr_1) delete(minute_hamsr_1) delete(hour_hamsr_2) delete(minute_hamsr_2) end do ; end loop over all 3 storms in 2016 ; average of all files avaps_temp_avg = dim_avg_n(avaps_temp_all,0) avaps_abshum_avg = dim_avg_n(avaps_abshum_all,0) hamsr_temp_avg = dim_avg_n(hamsr_temp_all,0) hamsr_abshum_avg = dim_avg_n(hamsr_abshum_all,0) hamsr_spechum_avg = dim_avg_n(hamsr_spechum_all,0) avaps_spechum_avg = dim_avg_n(avaps_spechum_all,0) ; difference of the average of the two diff_temp_avg = avaps_temp_avg - hamsr_temp_avg diff_abshum_avg = avaps_abshum_avg - hamsr_abshum_avg diff_spechum_avg = avaps_spechum_avg - hamsr_spechum_avg ; plot average profiles ; temp first ;************************************************ ; create plots, Temperature, then RH ;************************************************ wks = gsn_open_wks ("png","avaps_hamsr_temp_avg_2016") ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "AVAPS and HAMSR AVG 2016" 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@tmXBLabelStride = 2 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,avaps_temp_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 plot1 = gsn_csm_xy(wks,hamsr_temp_avg,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,avaps_temp_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) plot4 = gsn_csm_xy(wks,hamsr_temp_avg,lev_hamsr_1,res) overlay(plot0,plot3) overlay(plot0,plot4) draw(plot0) frame(wks) ;;;;;;;;;;;; ; absolute humidity ;;;;;;;;;;;;; wks = gsn_open_wks ("png","avaps_hamsr_abshum_avg_2016") res = True ; plot mods desired res@tiMainString = "AVAPS and HAMSR AVG 2016"; + 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@tmXBLabelStride = 2 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,avaps_abshum_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 plot1 = gsn_csm_xy(wks,hamsr_abshum_avg,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,res) ; overlay(plot0,plot2) ;end do res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot3 = gsn_csm_xy(wks,avaps_abshum_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) plot4 = gsn_csm_xy(wks,hamsr_abshum_avg,lev_hamsr_1,res) overlay(plot0,plot3) overlay(plot0,plot4) draw(plot0) frame(wks) delete(res) ;;;;;;;;;;;;;;;;;;;;; ;;; specific humidity ;;;;;;;;;;;;;;;;;;;;; wks = gsn_open_wks ("png","avaps_hamsr_spechum_avg_2016") res = True ; plot mods desired res@tiMainString = "AVAPS and HAMSR AVG 2016"; + 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.1 res@trXMaxF = 10.0 ;res@trYLog = True res@trYMinF = 0.1 res@trYMaxF = 1000.0 ;res@tmXBLabelStride = 2 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 = "Specific Humidity (g/kg)" res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,avaps_spechum_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 plot1 = gsn_csm_xy(wks,hamsr_spechum_avg,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,res) ; overlay(plot0,plot2) ;end do res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot3 = gsn_csm_xy(wks,avaps_spechum_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) plot4 = gsn_csm_xy(wks,hamsr_spechum_avg,lev_hamsr_1,res) overlay(plot0,plot3) overlay(plot0,plot4) draw(plot0) frame(wks) ;;;;;;;;;;;;;;;;;;;;; ; difference plots of the average ;;;;;;;;;;;;;;;;;;;;;; ; plot average profiles ; temp first ;************************************************ ; create plots, Temperature, then RH ;************************************************ wks = gsn_open_wks ("png","avaps_hamsr_temp_diff_2016") ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "Temperature" 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@gsnXRefLine = 0.0 res@gsnXRefLineColor = "Black" res@gsnXRefLineDashPattern = 2 res@gsnXRefLineThicknessF = 2 res@trXMinF = -2 res@trXMaxF = 2 ;res@trYLog = True res@trYMinF = 0.1 res@trYMaxF = 1000.0 ;res@tmXBLabelStride = 2 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 = "AVAPS - HAMSR (degC)" res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,diff_temp_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 draw(plot0) frame(wks) ;;;;;;;;;;;; ; absolute humidity ;;;;;;;;;;;;; wks = gsn_open_wks ("png","avaps_hamsr_abshum_diff_2016") ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "Absolute Humidity" 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 = -2 res@trXMaxF = 2 ;res@trYLog = True ;res@tmXBLabelStride = 2 res@trYMinF = 0.1 res@trYMaxF = 1000.0 res@gsnXRefLine = 0.0 res@gsnXRefLineColor = "Black" res@gsnXRefLineDashPattern = 2 res@gsnXRefLineThicknessF = 2 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 = "AVAPS - HAMSR (g/m3)" res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,diff_abshum_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 draw(plot0) frame(wks) ;;;;;;;;;;;;;;;;;;;;; ;;; specific humidity ;;;;;;;;;;;;;;;;;;;;; wks = gsn_open_wks ("png","avaps_hamsr_spechum_diff_2016") ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "Specific Humidity" 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@tmXBLabelStride = 2 res@trXMinF = -2 res@trXMaxF = 2 ;res@trYLog = True res@gsnXRefLine = 0.0 res@gsnXRefLineColor = "Black" res@gsnXRefLineDashPattern = 2 res@gsnXRefLineThicknessF = 2 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 = "AVAPS - HAMSR (g/kg)" res@xyLineThicknessF = 2.5 res@xyLineColor = colors_plot(0) plot0 = gsn_csm_xy(wks,diff_spechum_avg,lev_hamsr_1,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 1.0 draw(plot0) frame(wks) end