load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" begin ;******************** READ WRF ****************************************** ;WRf files all_files = systemfunc ("ls /data14a/maurom/SIMULATIONS/D_17/domain2/wrfout_d01_2017*") ;print(all_files) ;exit simulation = "D17_d02" longname = "WS_WD_10m" compound1 = "U10" compound2 = "V10" ; Remember to change the variables ncar_rea1 = "uwnd.10m.gauss.2017.nc" ncar_rea2 = "vwnd.10m.gauss.2017.nc" compound = "WIND" title = compound title1 = compound + "_" + simulation diro = title1 system("mkdir -p " + diro) ;************************************************************************ ; The WRF ARW input file. a = addfiles (all_files, "r") b = addfile ("/data14a/maurom/SIMULATIONS/D_17/domain2/wrfout_d01_2017-01-01_00:00:00","r") b1 = addfile(ncar_rea1,"r") b2 = addfile(ncar_rea2,"r") ; Check Variables ;print(a) ;print(b1) ;exit ;************************************************************************ ; Variables wr1 = wrf_user_getvar(a,compound1,-1) ; printVarSummary(wr1) wr2 = wrf_user_getvar(a,compound2,-1) ; printVarSummary(wr2) nc1 = b1->uwnd(:,:,:) nc1 = lonFlip(nc1) ; change longitude (0;360) -> (-180;180) ; printVarSummary(nc1) nc2 = b2->vwnd(:,:,:) nc2 = lonFlip(nc2) ; change longitude (0;360) -> (-180;180) ; printVarSummary(nc2) ; exit ;************************************************************************ ;Times Variable ; times = wrf_user_getvar(a,"times",-1) ; get all times in the file times = wrf_times_c(a[:]->Times,0) ntimes = dimsizes(times) ; number of times in the file ti = b1->time(:) dimtime = dimsizes(ti) ; number of times in the avarged var ntime = dimtime(0) ; Array to hold month abbreviations. Don't store anything in index ; '0' (i.e. let index 1=Jan, 2=Feb, ..., index 12=Dec). month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \ "Oct","Nov","Dec"/) ; Convert to UTC time. utc_date = cd_calendar(ti, 0) ; Store return information into more meaningful variables. year = tointeger(utc_date(:,0)) ; Convert to integer for month = tointeger(utc_date(:,1)) ; use sprinti day = tointeger(utc_date(:,2)) hour = tointeger(utc_date(:,3)) minute = tointeger(utc_date(:,4)) second = utc_date(:,5) ; Write out strings in the format "hhZ dd mmm yyyy". ; date_str = sprinti("%0.2iZ ", hour) + sprinti("%0.2i ", day) + \ date_str = sprinti("%0.2i ", day) + month_abbr(month) + " " + sprinti("%0.4i", year) ; print(date_str) ; exit ;************************************************************************ ; Daily Average ; wr&Time= wrf_times_c(a->Times,0) wr1&Time = times wr2&Time = times ; printVarSummary(wr_plane&Time) ; printMinMax(cd_calendar(wr_plane&Time,3),0) ; print(cd_calendar(wr_plane&Time,3)) ;************************************************************************** ; Variables to print wr1_avg = calculate_monthly_values(wr1, "avg", 0, False) wr1_avg@long_name= "WRFChem_U10" wr1_avg@units = "[m/s]" wr2_avg = calculate_monthly_values(wr2, "avg", 0, False) wr2_avg@long_name= "WRFChem_V10" wr2_avg@units = "[m/s]" delete(wr1) delete(wr2) nc1_avg = calculate_monthly_values(nc1, "avg", 0, False) nc1_avg@long_name = "NCAR/NCEP_U10" nc1_avg@units = "[m/s]" nc2_avg = calculate_monthly_values(nc2, "avg", 0, False) nc2_avg@long_name = "NCAR/NCEP_V10" nc2_avg@units = "[m/s]" delete(nc1) delete(nc2) wiwr_sp = sqrt (wr1_avg*wr1_avg+wr2_avg*wr2_avg) copy_VarCoords(wr1_avg, wiwr_sp) wiwr_sp@long_name= "WRFChem_Wsp" wiwr_sp@units = "[m/s]" ;delete(wr1_avg) ;delete(wr2_avg) winc_sp = sqrt (nc1_avg*nc1_avg+nc2_avg*nc2_avg) copy_VarCoords(nc1_avg, winc_sp) winc_sp@long_name= "NCAR/NCEP_Wsp" winc_sp@units = "[m/s]" ;delete(nc1_avg) ;delete(nc2_avg) ; Check Variables ; printVarSummary(wr1_avg) ; printVarSummary(wr2_avg) ; printVarSummary(nc1_avg) ; printVarSummary(nc2_avg) ; printVarSummary(wiwr_sp) ; printVarSummary(winc_sp) ;exit ;************************************************************************ ; We generate plots, but what kind do we prefer? ; type = "x11" type = "png" ; type = "pdf" ; type = "ps" ; type = "ncgm" type@wkWidth = 2500 type@wkHeight = 2500 ; wks = gsn_open_wks(type,"EU_d01") wks = gsn_open_wks(type,diro + "/" + compound) ;************************************************************************ ;Single plot options res = True res@gsnDraw = False ; Don't draw plots res@gsnFrame = False ; Don't advance frame res = wrf_map_resources(b,res) ; Add necessary resources for WRF map res@gsnAddCyclic = False ; Don't add longitude cyclic point res@cnFillOn = True ; Turn on color res@lbOrientation = "Horizontal" ; Labelbar ; res@lbLabelBarOn = False ; Turn off individual labelbars res@lbLabelBarOn = False ; Turn off individual labelbars res@cnLineLabelsOn = False ; label contour res@cnConstFLabelFontHeightF = 0.0 ; Costant field ; Map contour ; res@mpGeophysicalLineColor = "Black" ; res@mpGridLineColor = "Black" ; res@mpLimbLineColor = "Black" ; res@mpNationalLineColor = "Black" ; res@mpPerimLineColor = "Black" ; res@mpUSStateLineColor = "Black" ; mnmxint = nice_mnmxintvl( min((/min(wr_geo),min(nc_geo)/)), max((/max(wr_geo),max(nc_geo)/)), 18, False) res@cnLevelSelectionMode = "ManualLevels" ; Set contour levels to the same levels ; res@cnMinLevelValF = mnmxint(0) ; for both plots. ; res@cnMaxLevelValF = mnmxint(1) ; res@cnLevelSpacingF = mnmxint(2) res@cnMinLevelValF = 0 ; for both plots. res@cnMaxLevelValF = 8 res@cnLevelSpacingF = 1 res@cnFillPalette = "wind_17lev" ; res@cnFillPalette = "GMT_haxby" ; res@cnFillPalette = "sunshine_9lev" ; Set contour color ; res@cnLinesOn = True ; Contour Lines on/off res@cnLinesOn = False ; Contour Lines on/off ; res@cnLineColor = "Navy" ; res@mpPerimLineColor = "Black" ; res@mpPerimLineFontHeight = "Black" ; res@mpUSStateLineColor = "Black" res@mpGeophysicalLineColor = "Black" ; color of cont. outlines res@mpGeophysicalLineThicknessF = 1.5 ; thickness of outlines ;************************************************************************ ; Resources for vectors vecres = True ; vector only resources vecres@gsnDraw = False ; don't draw vecres@gsnFrame = False ; don't advance frame vecres@vcGlyphStyle = "LineArrow" ; curly vectors vecres@vcRefMagnitudeF = 5 ; define vector ref mag vecres@vcRefLengthF = 0.045 ; define length of vec ref vecres@gsnRightString = " " ; turn off right string vecres@gsnLeftString = " " ; turn off left string vecres@tiXAxisString = " " ; turn off axis label vecres@vcMinDistanceF = 0.03 ; turn off axis label ; vecres@vcRefAnnoOrthogonalPosF = -.535 ; move ref vector into plot ; vecres@vcMonoLineArrowColor = False ; color arrows based on magnitude ; vecres@vcRefLengthF = 0.03313608 vecres@vcMinFracLengthF = 0.3 ;********************************************************************** ;************************************************************************ ;---Set resources specific to WRF and NCEP plots res_wrf = res res_nc = res res_wrf@tfDoNDCOverlay = True ; Use the native WRF map projection ;********************************************************************** ; Find the time intersection between the WRF and NCEP arrays. ; - First convert WRF's time to be the sames units/calendar. ; - Use venn2_intersection to find where the two arrays have the same ; time values. ; - Use this new time array and coordinate subscripting to get the ; correct time index for each array for plotting. ;********************************************************************** ;---Convert WRF var's time to be the same as NCEP's var's time wiwr_sp&Time = cd_convert(wiwr_sp&Time,winc_sp&time@units) wr1_avg&Time = cd_convert(wr1_avg&Time,nc1_avg&time@units) wr2_avg&Time = cd_convert(wr2_avg&Time,nc2_avg&time@units) ;---Find the intersection between the two time arrays time_intersect = venn2_intersection(wiwr_sp&Time,winc_sp&time) copy_VarAtts(winc_sp&time,time_intersect) time_intersect_v1 = venn2_intersection(wr1_avg&Time,nc1_avg&time) copy_VarAtts(nc1_avg&time,time_intersect_v1) time_intersect_v2 = venn2_intersection(wr2_avg&Time,nc2_avg&time) copy_VarAtts(nc2_avg&time,time_intersect_v2) ; printVarSummary(time_intersect) ; printVarSummary(time_intersect_v1) ; printVarSummary(time_intersect_v2) ;exit ; Note: it's not necessary to subset nc and wr_avg and save ; to new variables like I'm doing here. I'm mainly doing this ; for illustrative purposes, so you can see what subsetting ; an array by coordinate values looks like. ; wr_avg_subset = wiwr_sp({time_intersect},:,:) nc_avg_subset = winc_sp({time_intersect},:,:) wr1_avg_subset = wr1_avg({time_intersect_v1},:,:) wr2_avg_subset = wr2_avg({time_intersect_v2},:,:) nc1_avg_subset = nc1_avg({time_intersect_v1},:,:) nc2_avg_subset = nc2_avg({time_intersect_v2},:,:) ; printVarSummary(wr_avg_subset) ; printVarSummary(nc_avg_subset) yyyymmddhh = cd_calendar(time_intersect,1) ; For nicer titles, if needed ntime_subset = dimsizes(time_intersect) ;print(yyyymmddhh) ;print(time_intersect) ;print(ntime_subset) ;exit ;********************************************************************** ; Loop over average time do it=0,ntime_subset-1 plot_wrg = gsn_csm_contour_map(wks,wr_avg_subset(it,:,:),res_wrf) plot_wrv = gsn_csm_vector(wks,wr1_avg_subset(it,:,:),wr2_avg_subset(it,:,:),vecres) overlay(plot_wrg,plot_wrv) ; result will be plotA plot_wrf = plot_wrg ; now assign plotA to array plot_ncg = gsn_csm_contour_map(wks,nc_avg_subset(it,:,:),res_nc) plot_ncv = gsn_csm_vector(wks,nc1_avg_subset(it,:,:),nc2_avg_subset(it,:,:),vecres) overlay(plot_ncg,plot_ncv) ; result will be plotA plot_ncr = plot_ncg ; now assign plotA to array ;---Create a panel of plots with 2 rows and 2 columns. pres = True ; Set panel resources. pres@gsnMaximize = True ; Maximize plots in panel. pres@pmLabelBarWidthF = 0.9 ; Change labelbar width. pres@gsnPanelLabelBar = True ; Turn on panel labelbar. ; pres@gsnPanelMainString = title + " : " + yyyymmddhh(it) pres@gsnPanelMainString = "Wind Speed and Wind direction "+yyyymmddhh(it) gsn_panel(wks,(/plot_wrf,plot_ncr/),(/1,2/),pres) end do ;---ImageMagick command cmd = "convert -delay 70. -loop 0. " + diro + "/*.png " + diro + "/" + title + ".gif" system(cmd) delete(compound) delete(wrfout) delete(wr) delete(nc) delete(times) delete(time_intersect) delete(wr_avg_subset) delete(nc_avg_subset) delete(yyyymmddhh) delete(ntime_subset) end