function main(arg) model='HRRR' param = 'simulated_radar_pnw' param2= '' map_save_loc = '/home/wxnorthw/public_html/weather/grads_maps/' bkg_map_loc = '/home/wxnorthw/bin/contents/resources/backgrounds/' if (arg='') say 'usage: run run file_name (without the .gs) YYYYMMDDHH ' say ' where ' say ' YYYY is the year' say ' MM is the month (01, 02, ... 12)' say ' DD is the day (01, 02, ... 31)' say ' HH is the hour(00, 06, 12 or 18) and ' * say ' LLLL is level in mb' say '' return endif *------ reinitialize --------------------------------------- say '' 'reinit' 'set xsize 800 600' *------ set the level ------------------------------------- say '' *level = subwrd(arg, 2) *say 'Level: ' level *------ set the model run date ---------------------------- mrundat = subwrd(arg, 1) mrundat = substr(mrundat, 1,8) say 'Model date: ' mrundat *------ set the model run hour ---------------------------- mruntim = subwrd (arg, 1) mruntim = substr(mruntim,9,2) say 'Model run time: ' mruntim 'Z' *------ get the model run day of month -------------------- mrundom = substr(mrundat, 7,2) *------ get the model run month number and name ----------- mrunmon = substr(mrundat, 5,2) mrunmn = num2mon(mrunmon) *------ get the model run year ---------------------------- mruny = substr(mrundat, 3,2) say '' *------ clear the display ---------------------------------- 'clear' *------ get info from nomad ------------------------------- 'sdfopen http://nomads.ncep.noaa.gov:9090/dods/hrrr/hrrr' mrundat '/hrrr_sfc_' mruntim 'z' *------ turn off frame, x/y labels ------------------------- 'set frame off' 'set xlab off' 'set ylab off' 'set parea 1 10 .3 6.7 ' 'set background 0' *------ looping message ------------------------------------ say '==========================================================' say ' Begin Processing Timesteps' say '==========================================================' say '' *------ loop through time steps ---------------------------- count = 0 ;while (count < 16); count = count+1 'set t ' count *------ get a formatted count ------------------------------ fcount = math_format('%02.0f', count) say 'Processing time step 'fcount *------ get forecast hour ---------------------------------- fcsthour = (count-1) fcsthour = math_format('%02.0f', fcsthour) *------ plot the display ----------------------------------- 'clear' ;*--- set special colors for this map... 'set rgb 17 255 255 255' 'set rgb 18 150 37 22';* map background... 'set rgb 19 0 43 215';* wind barbs... 'set rgb 21 255 255 0' 'set rgb 22 255 215 0' 'set rgb 23 205 133 0' 'set rgb 24 255 127 0' 'set rgb 25 238 64 0' 'set background 0' ;*--- set map coverage... 'set mproj scaled' 'set mpdset hires' 'set map 18 1 2' 'set lat 41.5 50' 'set lon -128 -113.8' 'set background 0' ;*--- turn off grads copyright and set map parameters... 'set grads off' ;*--- do radar... 'set gxout shaded' 'set rgb 39 4 233 231' 'set rgb 40 1 159 244' 'set rgb 41 3 0 244' 'set rgb 42 25 206 39' 'set rgb 43 0 142 0' 'set rgb 44 253 248 2' 'set rgb 45 255 148 0' 'set rgb 46 255 148 0' 'set rgb 47 253 0 0' 'set rgb 48 188 0 0' 'set rgb 49 188 0 0' 'set rgb 50 248 0 253' 'set rgb 51 152 98 198' 'set rgb 52 255 255 255' 'set clevs 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75' 'set ccols 1 39 40 41 42 43 44 45 46 47 48 49 50 51 52 ' 'display refcclm' *----------------------------------------------------------- * do the titles - * - *----------------------------------------------------------- 'set t ' count timestr = subwrd(result,4) vtz_ele = c2s(timestr) ;*--- get zulu time stuff... vh = subwrd(vtz_ele, 4) vh = math_format('%02.0f',vh) vd = subwrd(vtz_ele, 3) vd = math_format('%02.0f',vd) vm = subwrd(vtz_ele, 2) vm = math_format('%02.0f',vm) vmn = num2mon(vm) vy = subwrd(vtz_ele,1) vy = substr(vy,3,2) ;*--- get local time stuff... timestr = local_time(timestr) vtl_ele = c2s(timestr) vhl = subwrd(vtl_ele, 4) vhl = 1*vhl if(vhl>=12) ampm = 'pm' else ampm = 'am' endif if(vhl>12) vhl = vhl-12 endif if(vhl=0) vhl = '12' endif say 'hour: ' vhl vdl = subwrd(vtl_ele, 3) vml = subwrd(vtl_ele, 2) vmln = num2mon(vml) vyl = subwrd(vtl_ele, 1) vdln = subwrd(day_of_week(vyl ':' vml ':' vdl),2) ;*--- make color 99 dark gray for titles 'set rgb 99 64 64 64' 'set font 4' ;*------ local valid time... 'set string 99 bl 3' 'set strsiz .15' 'draw string 3.2 7.055 ' vhl ampm ' ' vdln ' ' vdl ' ' vmln ;*--- plot forecast hour... fh = fcsthour 'draw string 7.3 7.055 ' fcsthour ;*--- plot Z valid time and model run time... 'set string 15 bl 1' 'set strsiz .10' 'draw string 1.0 .08 (c) Weather NorthWest' 'draw string 4.2 .08 Valid Time: ' vh 'Z ' vd vmn vy 'draw string 7.5 .08 Model Run Time: 'mruntim 'Z ' mrundom mrunmn mruny *----------------------------------------------------------- * end title routine - *----------------------------------------------------------- *------ save the file -------------------------------------- say 'Saving file 'model'_'param'_'fcount'.png' say '' * 'printim ' map_save_loc model '_' param '_' fcount '.png x800 y600' 'printim ' map_save_loc model '_' param '_' fcount '.png x800 y600 -b -t 0' bkg_map_loc 'hrrr_nw_oregon_radar_background.png' endwhile return *----------------------------------------------------------- * * function c2s (colon to space) * * replaces colons with spaces in a string * *----------------------------------------------------------- function c2s(str) retstr='' count=0; while(count=9) gmt_offset = -8 endif endif if(day>1*first_sunday) gmt_offset = -8 endif endif if (month = 3) ;* march month_start_day = subwrd(day_of_week(year ':03:1:00'),3) second_sun_str = '8 14 13 12 11 10 9' second_sunday = subwrd(second_sun_str, month_start_day) if(1*day=1*second_sunday) if(hour<9) gmt_offset = -8 endif endif if(day<1*second_sunday) gmt_offset = -8 endif endif ;*--- calculate local time... local_hour = hour + gmt_offset local_day = day local_month = month local_year = year if(local_hour<0) local_hour = local_hour + 24 local_day = day - 1 if(local_day<1) local_month = month - 1 if(local_month<1) local_month = 12 local_year = year-1 endif if(local_month=1);local_day = 31; endif if(local_month=2) if(math_fmod(local_year,4)=0) local_day = 29 else local_day = 28 endif if(local_year=2000) local_day = 28 endif endif if(local_month=3);local_day = 31; endif if(local_month=4);local_day = 30; endif if(local_month=5);local_day = 31; endif if(local_month=6);local_day = 30; endif if(local_month=7);local_day = 31; endif if(local_month=8);local_day = 31; endif if(local_month=9);local_day = 30; endif if(local_month=10);local_day = 31; endif if(local_month=11);local_day = 30; endif if(local_month=12);local_day = 31; endif endif endif ;*--- set the timezone abbreviation... if(gmt_offset = -7) timezone='PDT' else timezone='PST' endif return local_year ' ' math_format('%02.0f', local_month) ' ' math_format('%02.0f', local_day) ' ' math_format('%02.0f', local_hour) ' ' timezone *----------------------------------------------------------- * * function day_of_week(grads_time) * * returns 0-6 (Sun thru Sat) given a date in the format * yyyy:mm:dd:hh (or yyyy:mm:dd - the :hh is ignored) * * uses Zellers algorithm from wikipedia * * returns: day_name day_abbreviation day_of week * (ie Sunday Sun 1) * where day_of_week is 1 for Sun, 2 for Mon, etc. * *----------------------------------------------------------- function day_of_week(time_string) ;*--- set the day names... day_abbr_string = "Sat Sun Mon Tue Wed Thu Fri" day_name_string = "Saturday Sunday Monday Tuesday Wednesday Thursday Friday" day_num_string = "7 1 2 3 4 5 6" ;*--- get variables needed for calculations... time_string = c2s(time_string) d = subwrd(time_string, 3) m = subwrd(time_string, 2) Y = subwrd(time_string, 1) if(m<3) Y=Y-1 m = m+12 endif y = substr(Y, 3, 2) c = substr(Y, 1, 2) ;*--- do the calculations and return a day name... w = math_fmod((d + math_int((13*(m + 1))/5) + y +math_int(y/4) + math_int(c/4) + 5*c),7)+1 day_abbr = subwrd(day_abbr_string, w) day_name = subwrd(day_name_string, w) day_num = subwrd(day_num_string, w) return day_name ' ' day_abbr ' ' day_num *------ script end -----------------------------------------