load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/calendar_decode2.ncl" begin ; 1. Read file and Print Meta Info and File name is being passed on to the code using the Shell script ; name_of_file = "21apr2010.ua.nc" SPLIT_FILE = str_split(name_of_file,".") print(SPLIT_FILE) ; storms_file = "dcc.wcc."+SPLIT_FILE(0)+".txt" storms_file = "dcc."+SPLIT_FILE(0)+".txt" wcc40_file = "wcc."+SPLIT_FILE(0)+".txt" ; storms_file = "dcc.wcc.2008-04-20.txt" print(storms_file) dcc_which_are_also_wcc = "dcc.wcc."+SPLIT_FILE(0)+".txt" ;;;; ; LOCATIONS of nearby RS stations File = addfile(name_of_file,"r") var_names = getfilevarnames(File) print(var_names) print(File) lon = File->longitude lat = File->latitude lon@units = "degrees_east" lat@units = "degrees_north" ; printVarSummary(lon) ; printVarSummary(lat) Lev = File->level print(Lev) Dumy = File->time ; print(Dumy) Time = calendar_decode2(Dumy,3) ; if 3, O/P 20110414, if 2 O/P 2011041400 and so on ; printVarSummary(Time) NVars = getfiledimsizes(File) print(NVars) Var_names = getvardims(File) Ntimes = dimsizes(Time) ; print(Ntimes+" "+Var_names+" "+Ntimes) ; Extracting variables from File U_shrt = File->u U_flt = short2flt(U_shrt) ; printVarSummary(U_flt) u_flt = U_flt(0,5,:,:) V_shrt = File->v V_flt = short2flt(V_shrt) ; printVarSummary(V_flt) v_flt = U_flt(0,5,:,:) Q_shrt = File->q Q_flt = 1000*short2flt(Q_shrt) ; 1000 is multiplied for converting to g/kg ; Q_flt = smth9_Wrap(Q_FLT,0.5,0.25,False) ; printVarSummary(Q_flt) ; printMinMax(Q_flt,True) ;;; Extracting variables from file is OVRRRRRRRRRRRRRRRR base = new(4,graphic) plot = new(4,graphic) plot_q = new(4,graphic) add_DCC_WCC = new(4,graphic) add_DCC40_WCC40 = new(4,graphic) add_RS_STN_LOC = new(4,graphic) add_WCC40 = new(4,graphic) india_id = new(4,graphic) ; plot1 = new(4,graphic) res = True ; Set -----------MAP ---- res@mpFillOn = False res@mpLimitMode = "Corners" res@mpLeftCornerLatF = 5.0 res@mpLeftCornerLonF = 60.0 res@mpRightCornerLatF = 30.0 res@mpRightCornerLonF = 100.0 ; ::::::: TICK MARK:::: res@tmXTOn = False res@tmYROn = False res@tmXBMode = "Explicit" ; res@tmXBValues = (/"60","65","70","75","80","85","90","95","100"/) ; res@tmXBLabels = (/"60~S~o~N~E","65~S~o~N~E","70~S~o~N~E","75~S~o~N~E","80~S~o~N~N","85~S~o~N~E","90~S~o~N~E","95~S~o~N~E","100~S~o~N~E"/) res@tmXBMajorLengthF = 0.0 res@tmXBMinorLengthF = 0.0 res@tmXBMinorOn = False res@tmYLMode = "Explicit" ; res@tmYLValues = (/"5","10","15","20","25","30"/) ; res@tmYLLabels = (/"5~S~o~N~N","10~S~o~N~N","15~S~o~N~N","20~S~o~N~N","25~S~o~N~N","30~S~o~N~N"/) res@tmYLMajorLengthF = 0.0 res@tmYLMinorLengthF = 0.0 res@tmYLMinorOn = False ;:::::: GSN resources:::::::::::::: res@gsnDraw = False res@gsnFrame = False res@gsnMaximize = True ;::::: Countour resources:::::::: resq = res delete(resq@mpFillOn) delete(resq@mpLimitMode) delete(resq@mpLeftCornerLatF) delete(resq@mpLeftCornerLonF) delete(resq@mpRightCornerLatF) delete(resq@mpRightCornerLonF) resq@cnFillOn = False resq@cnMonoFillColor = False ; resq@cnFillColors = False ; resq@cnFillPalette = "cbr_wet" resq@cnLevelSelectionMode = "ExplicitLevels" resq@cnLevels = ispan(4,20,4) resq@cnLineLabelsOn = True resq@cnLineLabelPlacementMode = "computed" resq@cnLineLabelDensityF = 2.0 ; will put more num of labels on the contour resq@cnLineLabelInterval = 1.0 ; will put labels on every line resq@cnInfoLabelOn = False ; will put off the rectangular box that will show the contour range resq@cnLinesOn = True resq@cnLineColor = "purple" resq@cnLineThicknessF = 2.5 resq@cnSmoothingOn = True resq@cnSmoothingTensionF = 1.5 ; will smooth the contour resq@cnLineLabelBackgroundColor = "pink" resq@cnLineLabelFontColor = "blue" res@cnFillOn = True res@cnMonoFillColor = False ; res@cnFillPalette = "WhBlGrYeRe" COLORS = (/"white","darkseagreen1","darkseagreen2","darkseagreen3","darkseagreen4"/) res@cnFillPalette = COLORS ; res@cnFillPalette = "GMT_topo" res@cnLevelSelectionMode = "ExplicitLevels" ; res@cnLevels = (/20,25,30/) ;(/10,15,20,30/); (/4,6,8,10/) depending on the selection of pressure level res@cnFillOn = True res@cnLineLabelsOn = False res@cnLinesOn = False ; delete(res@cnFillPalette) ; ::::::: LABEL Resources:::::::::::::::: res@lbLabelBarOn = False res@lbOrientation = "vertical" ; ---Adding Shape file for Indian Region---- IND_shape_fname = "IND_adm1.shp" shpres = True shpres@gsLineColor = "Black" shpres@gsLineThicknessF = 0.5 delim = " " RES_MARK_DCC = True RES_MARK_DCC@gsMarkerIndex = 13 RES_MARK_DCC@gsMarkerSizeF = 0.01 RES_MARK_DCC@gsMarkerColor = "PINK" ;DCC RES_MARK_DCC@gsMarkerThicknessF = 4.0 RES_MARK_WCC = RES_MARK_DCC RES_MARK_WCC@gsMarkerIndex = 14 RES_MARK_WCC@gsMarkerColor = "NavyBlue" ;" RES_MARK_DCC_WCC = RES_MARK_DCC RES_MARK_DCC_WCC@gsMarkerColor = "steelblue1" RES_MARK_DCC40_WCC40 = RES_MARK_DCC RES_MARK_DCC40_WCC40@gsMarkerColor = "red" RES_MARK_WCC40 = RES_MARK_DCC ; MARKING WCC40, color RES_MARK_WCC40@gsMarkerIndex = 14 RES_MARK_WCC40@gsMarkerColor = "palevioletred4" RES_RS = RES_MARK_DCC RES_RS@gsMarkerIndex = 1 RES_RS@gsMarkerColor = "Blue" ; Set Vector wind res_vc = True res_vc@gsnDraw = False res_vc@gsnFrame = False res_vc@vcGlyphStyle = "CurlyVector" res_vc@vcLineArrowColor = "red" res_vc@vcMinDistanceF = 0.020 ;res_vc@vcRefMagnitudeF = 20.0 ; 10 or 20 depending on preesure level res_vc@vcRefLengthF = 0.040 ;res_vc@vcMonoLineArrowColor = False res_vc@vcRefAnnoOrthogonalPosF = -1.0 ; will move the Ref vector to top Right res_vc@vcRefAnnoArrowLineColor = "black" res_vc@vcMinDistanceF = 0.04 ; will put some distance between vectors ; res@vcLineArrowHeadMinSizeF = 0.01 ; res@vcLineArrowHeadMaxSizeF = 0.01 ; res@vcFillArrowsOn = True ; res@vcFillArrowEdgeColor = "black" ; res@vcFillArrowFillColor = "green" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; level (0) 200 (1) 300 (2) 500 (3) 700 (4) 850 (5) 925 (6) 1000 do lev = 0, 6 wks = gsn_open_wks("png","ecmf."+name_of_file+"."+lev) ; wks = gsn_open_wks("png","ecmf."+lev) if ((lev.eq.0).or.(lev.eq.1).or.(lev.eq.2)) then res_vc@vcRefMagnitudeF = 20.0 ; 10 or 20 depending on preesure level res@cnLevels = (/20,25,30/) else delete(res@cnLevels) res_vc@vcRefMagnitudeF = 10.0 ; 10 or 20 depending on preesure level res@cnLevels = (/10,20,30/) end if do deltaT = 0,Ntimes-1 ; if (deltaT.eq.0) then ; res@tmXBValues = (/"50","120"/) ; res@tmXBLabels = (/"50~S~o~N~E","120~S~o~N~E"/) ; res@tmYLValues = (/"5","10","15","20","25","30"/) ; res@tmYLLabels = (/"5~S~o~N~N","10~S~o~N~N","15~S~o~N~N","20~S~o~N~N","25~S~o~N~N","30~S~o~N~N"/) ; else if (deltaT.eq.1) then ; delete(res@tmXBValues) ; delete(res@tmXBLabels) ; res@tmXBValues = (/"50","120"/) ; res@tmXBLabels = (/"50~S~o~N~E","120~S~o~N~E"/) ; delete(res@tmYLValues) ; delete(res@tmYLLabels) ; res@tmYLValues = (/"0","40"/) ; res@tmYLLabels = (/"0~S~o~N~E","40~S~o~N~E"/) ; res@tmYLOn = False ; else if (deltaT.eq.2) then ; res@tmXBOn = True ; res@tmYLOn = True ; else if (deltaT.eq.3) then ; res@tmYLOn = False ; res@tmXBOn = True ; end if ; end if ; end if ; end if wndspd = (/sqrt(U_flt(deltaT,lev,:,:)^2+V_flt(deltaT,lev,:,:)^2)/) wndspd!0 = "lat" wndspd!1 = "lon" wndspd&lon = lon wndspd&lat = lat wndspd@lat = "degrees_north" wndspd@lon = "degrees_east" wndspd@_FillValue = -32767 ; printVarSummary(wndspd) ; WNDSPD = (/sqrt(u_flt^2+v_flt^2)/) ; WNDSPD@_FillValue = -32767 ; WNDSPD!0 = "lat" ; WNDSPD!1 = "lon" ; WNDSPD&lon = lon ; WNDSPD&lat = lat ; WNDSPD@lat = "degrees_north" ; WNDSPD@lon = "degrees_east" Q_flt!0 = "Time" Q_flt!1 = "Lev" Q_flt!2 = "lat" Q_flt!3 = "lon" Q_flt&lon = lon Q_flt&lat = lat Q_flt@_FillValue = -32767 ; printVarSummary(Q_flt) ; print(deltaT+" "+Lev(lev)) ; printMinMax(Q_flt(deltaT,lev,:,:),True) ; ---------------Plot the wind speed Filled Contours and Specific humidity res@gsnLeftString = Time(deltaT) res@gsnRightString = Lev(lev)+" hPa" ; res@vpHeightF = 0.555583 ; res@vpWidthF = 0.888933 base(deltaT) = gsn_csm_contour_map(wks,wndspd,res) plot_q(deltaT) = gsn_csm_contour(wks,Q_flt(deltaT,lev,:,:),resq) india_id(deltaT) = gsn_add_shapefile_polylines(wks,base(deltaT),IND_shape_fname,shpres) overlay(base(deltaT),plot_q(deltaT)) DCC_WCC = storms_file print(DCC_WCC) TSTR_DCC_WCC = asciiread(DCC_WCC,-1,"string") print(TSTR_DCC_WCC) LATC_DCC_WCC = stringtofloat(str_get_field(TSTR_DCC_WCC,1,delim)) LATC_DCC_WCC@_FillValue = -99.99 LONC_DCC_WCC = stringtofloat(str_get_field(TSTR_DCC_WCC,2,delim)) LONC_DCC_WCC@_FillValue = -99.99 add_DCC_WCC(deltaT) = gsn_add_polymarker(wks,base(deltaT),LONC_DCC_WCC,LATC_DCC_WCC,RES_MARK_DCC_WCC) ;;;;; DCC40 which are alsoWCC40 DCC40_which_are_also_WCC40 = dcc_which_are_also_wcc TSTR_DCC40_WCC40 = asciiread(DCC40_which_are_also_WCC40,-1,"string") print(TSTR_DCC40_WCC40) LATC_DCC40_WCC40 = stringtofloat(str_get_field(TSTR_DCC40_WCC40,1,delim)) LONC_DCC40_WCC40 = stringtofloat(str_get_field(TSTR_DCC40_WCC40,2,delim)) add_DCC40_WCC40(deltaT) = gsn_add_polymarker(wks,base(deltaT),LONC_DCC40_WCC40,LATC_DCC40_WCC40,RES_MARK_DCC40_WCC40) ;;;; WCC40 = wcc40_file print(WCC40) TSTR_WCC40 = asciiread(WCC40,-1,"string") print(TSTR_WCC40) LATC_WCC40 = stringtofloat(str_get_field(TSTR_WCC40,1,delim)) LATC_WCC40@_FillValue = -99.99 LONC_WCC40 = stringtofloat(str_get_field(TSTR_WCC40,2,delim)) LONC_WCC40@_FillValue = -99.99 add_WCC40(deltaT) = gsn_add_polymarker(wks,base(deltaT),LONC_WCC40,LATC_WCC40,RES_MARK_WCC40) ;;; Adding RS station locations RS_STN_LOC = "RS_stn_locations.txt" Readfile = asciiread(RS_STN_LOC,-1,"string") LAT_RS_STN = stringtofloat(str_get_field(Readfile,1," ")) LON_RS_STN = stringtofloat(str_get_field(Readfile,2," ")) add_RS_STN_LOC(deltaT) = gsn_add_polymarker(wks,base(deltaT),LON_RS_STN,LAT_RS_STN,RES_RS) ;:::: Vector wind :::::::: res_vc@gsnLeftString = "" res_vc@gsnRightString = "" ; res_vc@vpHeightF = 0.555583 ; res_vc@vpWidthF = 0.888933 plot(deltaT) = gsn_csm_vector(wks,U_flt(deltaT,lev,:,:),V_flt(deltaT,lev,:,:),res_vc) overlay(base(deltaT),plot(deltaT)) draw(base) ; delete(resq@cnLevels) ; getvalues base(0) ; get plot size for use ino ; "vpHeightF" : vph ; "vpWidthF" : vpw ; end getvalues ; print(vph+" base0 "+vpw) ; setvalues base(1) ; "vpHeightF" : vph ; "vpWidthF" : vpw ; end setvalues ; setvalues base(1) ; "vpHeightF" : vph ; "vpWidthF" : vpw ; end setvalues ; setvalues base(2) ; "vpHeightF" : vph ; "vpWidthF" : vpw ; end setvalues ; setvalues base(3) ; "vpHeightF" : vph ; "vpWidthF" : vpw ; end setvalues frame(wks) end do resP = True resP@gsnPanelLabelBar = True resP@lbOrientation = "Vertical" gsn_panel(wks,base,(/2,2/),resP) end do ; lev end