load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl" begin wks_type = "png" url = "http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg/" fnameprefix = "GFS_Global_0p5deg_" ; date_init = "12" fnamesuffix = "00.grib2/GC" level = 850 * 100 latmax = 45 latmin = -10 lonmax = 360 lonmin = 0 date_year = systemfunc("date +%Y") date_month = systemfunc("date +%m") date_day = "11"; systemfunc("date +%d") date_dayweek= systemfunc("date +%w") date_hour_s = systemfunc("date +%H") date_Zfactor= systemfunc("date +%:::z") print(date_hour_s) if (date_hour_s.eq."04") then date_init = "00" else if (date_hour_s.eq."10") then date_init = "06" else if (date_hour_s.eq."16") then date_init = "12" else if (date_hour_s.eq."22") then date_init = "18" else print("Ni pedo, usaremos esta:") date_init = "12" end if end if end if end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Get forecast times... ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; date_yearfull = systemfunc("date +%Y") yearint = tointeger(date_yearfull) monthint = tointeger(date_month) dayint = tointeger(date_day) daywint = tointeger(date_dayweek) dateint = tointeger(date_init) tStrt = ut_inv_calendar(yearint,monthint,dayint,dateint,0,0,"hours ref 1-1-1 00:00:0.0", 0) tstep = 3 ntim = 50 timeF = tStrt + ispan(0,tstep*(ntim-1),3) timeF@units = "hours since 1-1-1 00:00:0.0" ymdhF = cd_calendar(timeF, -3) ;printVarSummary(ymdhF) fchain = url + fnameprefix + date_year +date_month + date_day + "_" + date_init + fnamesuffix print(fchain) exists = isfilepresent(fchain) ; print(exists) if (exists.eq.True) then print("The file for time " + date_init + " exists") else print("The file for time " + date_init + " does not exist") end if exist2 = new(1,"logical") exist2 = exists ;print(exist2) vcounts = 1 do while((exist2.eq.False).and.(vcounts.lt.11)) sleep(180) exist2 = isfilepresent(fchain) print(exist2) vcounts = vcounts + 1 end do vdateend = systemfunc("date") print("File ready at " + vdateend) f = addfile (fchain, "r") vtime0 = f->time vtime1 = f->time1 vtime2 = f->time2 ; vtime3 = f->time3 print(dimsizes(vtime0)) print(dimsizes(vtime1)) print(dimsizes(vtime2)) ;print(dimsizes(vtime3)) if (vtime0(0).ne.0) then print("Salio en cero") vtime = doubletoint(vtime0) else if (vtime1(0).ne.0) then print("Salio en uno") vtime = doubletoint(vtime1) else if (vtime2(0).ne.0) then print("Salio en dos") vtime = doubletoint(vtime2) ; else if (vtime3(0).ne.0) then ; print("Salio en tres") ; vtime = vtime3 ; end if else print("There is nothing") vtime = doubletoint(vtime0) end if end if end if print(vtime0(0)) ;;;sleep(60) lat = f->lat lon = f->lon nlats = dimsizes(lat) nlons = dimsizes(lon) VU = f->$"u-component_of_wind_isobaric"$({vtime(0):vtime(49)},{level},{latmin:latmax},{lonmin:lonmax}) VU2 = VU / 0.5144 copy_VarMeta(VU,VU2) VV = f->$"v-component_of_wind_isobaric"$({vtime(0):vtime(49)},{level},{latmin:latmax},{lonmin:lonmax}) VV2 = VV / 0.5144 copy_VarMeta(VV,VV2) VT = f->Temperature_isobaric({vtime(0):vtime(49)},{level},{latmin:latmax},{lonmin:lonmax}) RH = f->Relative_humidity_isobaric({vtime(0):vtime(49)},{level},{latmin:latmax},{lonmin:lonmax}) MR2 = RH * 0.01 * 6.11 * (10^((7.5*(VT-273.15))/((VT-273.15)+237.3))) copy_VarMeta(VT,MR2) VG = f->Geopotential_height_isobaric({vtime(0):vtime(49)},{level},{latmin:latmax},{lonmin:lonmax}) VG2 = VG/10 copy_VarMeta(VG,VG2) VT2 = VT-273.15 copy_VarMeta(VT,VT2) ;*********************************************** ; Start the party... print("Starting the party...") do k = 0,49 print(k) vk = k + 2 wks_name = "/ct11/vtorres/webs/page/850/map_850_" + tostring(vk) vu2 = VU2(k,:,:) vv2 = VV2(k,:,:) vmr2 = MR2(k,:,:) vt2 = VT2(k,:,:) vg2 = VG2(k,:,:) ;printVarSummary(vmr2) print("Variables read...") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Map, title and size ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; wks_type@wkWidth = 1280 wks_type@wkHeight = 1024; 1500 wks = gsn_open_wks (wks_type, wks_name) ; open workstation ; gsn_define_colormap (wks,"CBR_wet") ; choose color map gsn_define_colormap (wks,"MPL_GnBu_MODvtp3") ; choose color map ; gsn_define_colormap (wks,"MPL_StepSeq") ; choose color map ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; All the resources here... ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; res = True ; plot mods desired vcres = True res@gsnMaximize = True ; res@gsnAddCyclic = False ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Contour resources ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; res@gsnDraw = False res@gsnFrame = False res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels res@cnMinLevelValF = 6. ; set min contour level res@cnMaxLevelValF = 20. ; set max contour level res@cnLevelSpacingF = 0.1; 1. ; set contour spacing res@lbLabelStride = 10; 2 ; skip every other label res@cnFillOn = True ; color fill res@cnLinesOn = False ; no contour lines res@gsnSpreadColors = True ; use total colormap res@gsnLeftString = "" res@gsnRightString = "" ; res@mpShapeMode = "FreeAspect" ; res@vpHeightF = 0.5 ; Changes the aspect ratio ; res@vpWidthF = 0.90 ; res@vpXF = 0.0 ; change start locations ; res@vpYF = 0.70 ; the plot ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; resT = True resT@gsnDraw = False resT@gsnFrame = False resT@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels resT@cnMinLevelValF = -21. ; set min contour level resT@cnMaxLevelValF = 21. ; set max contour level resT@cnLevelSpacingF = 3. ; set contour spacing resT@cnFillOn = False ; color fill resT@cnLinesOn = True ; no contour lines resT@cnLineThicknessF = 2.0 ; thicker lines resT@cnLineLabelsOn = True ; line labels resT@cnInfoLabelOn = False resT@cnLineLabelBackgroundColor = "Transparent" ; white bckgrnd around label resT@cnLineLabelAngleF = 0 resT@cnLineLabelFontHeightF = 0.008 ; resT@cnMonoLineLabelFontColor = False resT@cnLineLabelFontColor = "mediumvioletred" resT@gsnLeftString = "" resT@gsnRightString = "" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; resG = True resG@gsnDraw = False resG@gsnFrame = False resG@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels resG@cnMinLevelValF = 126. ; set min contour level resG@cnMaxLevelValF = 162. ; set max contour level resG@cnLevelSpacingF = 3. ; set contour spacing resG@cnLineLabelInterval = 1 resG@cnFillOn = False ; color fill resG@cnLinesOn = True ; no contour lines resG@cnMonoLineColor = True resG@cnLineColor = "midnightblue" resG@cnLineThicknessF = 3.0 resG@cnLineLabelsOn = True ; no line labels resG@cnInfoLabelOn = False resG@cnLineLabelFontColor = "midnightblue" resG@cnLineLabelBackgroundColor = "Transparent" ; white bckgrnd around label resG@cnLineLabelAngleF = 0 resG@cnLineLabelFontHeightF = 0.01 resG@gsnLeftString = "" resG@gsnRightString = "" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Label resources ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; res@lbBoxLinesOn = False res@lbTitleOn = True res@lbTitleString = "850hPa Mixing Ratio, Height, Temperature and Winds" res@lbTitleFontHeightF = 0.013 res@gsnLeftString = "Run: " + systemfunc("date +%a") + " " + systemfunc("date +%d-%b-%Y") + " " + date_init + "UTC";; res@gsnCenterString = "Forecast: " + vtime(k) + " hr" dateFv = ymdhF(k) ;print(dateFv) yyyy = dateFv/1000000 mmddhh = dateFv-yyyy*1000000 ; mmdd = yyyymmdd%10000 mm = mmddhh/10000 ddhh = mmddhh-mm*10000 ; dd = mmdd%100 dd = ddhh/100 hh = ddhh-dd*100 id_day= day_of_week(yyyy,mm,dd) month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", "Oct","Nov","Dec"/) day_tag = (/"","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18","19","20", \ "21","22","23","24","25","26","27","28","29","30","31"/) hour_tag = (/"03","","","06","","","09","","","12","","","15","","","18","","","21","","","00"/) day_nme_tag= (/"Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"/) ;print(yyyy) ;print(mm) ;print(dd) ;print(hh) ; res@gsnRightString = tostring(dateFv) + " UTC" res@gsnRightString = day_nme_tag(id_day) + " " + day_tag(dd) + "-" + month_abbr(mm) + "-" + tostring(yyyy) + " " + hour_tag(hh) + "UTC" res@gsnRightStringFontHeightF = 0.013 res@gsnCenterStringFontHeightF = 0.013 res@gsnLeftStringFontHeightF = 0.013 ;;; if (k.eq.0) then ;;; res@gsnRightStringFontColor = "black" ;;; res@gsnCenterString = "Initial time" ;;; else res@gsnRightStringFontColor = "red" res@gsnCenterStringFontColor = "red" ;;; end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Vector resources ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; vcres@gsnDraw = False vcres@gsnFrame = False vcres@gsnLeftString = " " vcres@gsnRightString = " " vcres@vcRefAnnoOn = False; True ; False if you dont want Reference vector vcres@vcRefMagnitudeF = 10; 20.0; 1.25 ; define vector ref mag vcres@vcRefLengthF = 0.018; 0.025; 0.045 ; define length of vec ref vcres@vcGlyphStyle = "WindBarb" ; turn on curly vectors vcres@vcMinDistanceF = 0.015; 0.01 ; vcres@vcLineArrowThicknessF = 1.0;1.5 ; vcres@vcFillArrowEdgeThicknessF = 1.0;2 ; vcres@vcLineArrowThicknessF = 1.0 ;1.5 vcres@vcRefAnnoString2 = "" vcres@vcLineArrowColor = "black" vcres@vcWindBarbColor = "gray30" vcres@vcRefAnnoString1 = vcres@vcRefMagnitudeF + " ms~S~-1~N~" vcres@vcRefAnnoOrthogonalPosF = -0.451; -1.40;-0.68 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Map resources ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; res@mpGridAndLimbOn = True ; turn on grid lines res@mpGridLatSpacingF = 10 res@mpGeophysicalLineThicknessF = 3 res@mpMinLatF = latmin res@mpMaxLatF = latmax res@mpMinLonF = -150 res@mpMaxLonF = -30 res@mpFillOn = False res@mpGeophysicalLineColor = "gray20" ; color of cont. outlines res@mpOutlineBoundarySets = "National"; turn on states res@mpDataBaseVersion = "mediumres" ; select database res@mpDataSetName = "Earth..2" res@pmTickMarkDisplayMode = "Always" ; turn on fancy tickmarks ; res@pmLabelBarWidthF = 0.4 ; control size of colorbar res@pmLabelBarHeightF = 0.075 ; res@pmLabelBarOrthogonalPosF = 0.15 ; position wrt plot res@lbLabelFontHeightF =.009 ; make labels larger ;res@tiMainString = "Western Region" res@tmXBLabelFontHeightF = 0.010 ; change maj lat tm spacing ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Draw map ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; plot = gsn_csm_contour_map(wks,vmr2,res) getvalues plot "mpAreaNames" : area_names end getvalues res@mpOutlineSpecifiers = area_names(617:655) plot = gsn_csm_contour_map(wks,vmr2,res) plot_vs = gsn_csm_vector(wks,vu2,vv2,vcres) plotT = gsn_csm_contour(wks,vt2,resT) plotT = ColorNegDashZeroPosContour(plotT,"blue","transparent","mediumvioletred") plotG = gsn_csm_contour(wks,vg2,resG) overlay(plot,plot_vs) overlay(plot,plotT) overlay(plot,plotG) draw (plot) frame(wks) dir_clocks = "/ct11/vtorres/webs/page/clocks3/" if (wks_type.eq."png") then system("convert -trim -density 400x400 "+ wks_name + ".png " + wks_name + ".png") system("composite -geometry +934+542 " + dir_clocks + "clock_" + hour_tag(hh) + "Z.png " + wks_name + ".png " + wks_name + ".png") else print("Ya estas peinado p'atras") end if ;; if (k.eq.0) then ;; system("cp " + wks_name + ".png " + wks_name + "_icon.jpg") ;; system("convert -resize 85% " + wks_name + "_icon.jpg " + wks_name + "_icon.jpg") ;; system("mv " + wks_name + "_icon.jpg /ct11/vtorres/webs/page/icons/") ;; else ;; end if end do print("Done it") print("***** Done 850 level ******") end