; Script to extract data from a coarse grid to plot for the nest(s) ; We want to extract inforomation just within the coordinates defined for the nests only. ; NOT for the whole domain...! load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;load "./contour_resources.ncl" load "./WRFCHEM_misc.ncl" begin DATADir = "ls -1 /Volumes/Seagate\ Bac/Full_wrfchem_run_frm_26" savepath = "grid_data_plots/" ;"../plots/height_interp_traj/" ;ncodes = dimsizes(runcode) ; ; domain = "01" ; parent grid (domain 1) ; domain = "02" ; 12t nest (domain 2 ; domain = "03" ; 2nd nest (domain 3) ; domain = "04" ; 3rd nest (domain 4) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Grid_list = (/"d01","d02","d03","d04"/) grid_list = (/"d02"/) ;ngrid = dimsizes(grid_list) ; define data directory DATAPath = ("../"+DATADir+"/") ; datapath = ("ls -1 /Volumes/Segate\ Bac/Full_wrfchem_run_frm_26") FILES = systemfunc (" ls -1 " + DATAPath + "wrfout_d01_2008-08-*") numFILES = dimsizes(FILES) specieslist = (/"o3"/) ; nspecies = dimsizes(specieslist) ;;;;;;loop through species for plotting do ispec = 0, nspecies-1 ; species species = specieslist(ispec) print("Working on " +species) ;;Define save name, path and open workstation. name = species+"_conctr'dx" path = "../trial_plots/concentrations/" wks = gsn_open_wks("pdf",path+name) gsn_define_colormap(wks,"WhiteBlueGreenYellowRed") ;------------------------------------------------------------------------------------- ;; Set some basic resources res = True mpres = True ; Map resources opts = True mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpLimitMode = "LatLon" mpres@mpGeophysicalLineThicknessF = 3.5 mpres@tmXBLabelFontHeightF = 0.014 mpres@tmYLLabelFontHeightF = 0.014 ;res@InitTime = False ;res@Footer = False ;res@CommonTitle = False res@MainTitle = "Concentration" res@cnFillOn = True res@cnLinesOn = True ;opts@cnFillMode = "AreaFill" opts@pmLabelBarOrthogonalPosF = -0.07 opts@cnMissingValFillColor = 20 opts@cnMissingValFillPattern = 11 pltres = True ; Plot resources pltres@FontHeightF = 0.02 ; set the font size of the titles pltres@CommonTitle = True ;pltres@FramePlot = False ; if set, will plot data in 1 graph for all time files res@trXMinF = 0.0 lnres = True lnres@gsLineThicknessF = 15.0 ; line size lnres_back =True ; Background line lnres_back@gsLineThicknessF = 18.0 lnres_back@gsLineColor = "black" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; firstime = True dt = 6 ;6 hr output do ifil = 0, numFILES-1, dt ; file loop a = addfile(FILES(ifil)+".nc","r") ; Open the next file ;................................................................................................... ;;What times and how many time steps are in the data set? times = wrf_user_list_times(a) ; get times in the file ntimes = dimsizes(times) ; number of times in the file xlat = wrf_user_getvar(a,"XLAT",0) xlon = wrf_user_getvar(a,"XLONG",0) ;ter = wrf_user_getvar(a, "HGT",0) do it = 0,ntimes-1 ; time loop ;print("Working on time: " + times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots ;;Specify the coordinate points (grid location) to plot ; d02, d03, d04 are subregions within d01. Define the coord. (corner) points for each grid on the ; parent domain ;;Define coord for d02 Max_Lat = -3.8 Min_Lat = -11.5 Max_Lon = 17.5 Min_Lon = 9.5 ;;Define coord for d03 ;Max_Lat = -3.8 ;Min_Lat = -14.0 ;Max_Lon = 44.5 ;Min_Lon = 33.0 ;;Define coord for d04 ;Max_Lat = -17.0 ;Min_Lat = -28.0 ;Max_Lon = 40.0 ;Min_Lon = 27.0 ;delete(griddata) ;llmin = (/Min_Lon,Min_Lat/) ; send info to arrays ;llmax = (/Max_Lon,Max_Lat/) ;;;;; Extract data for grid cells print("Extracting grid data for grid: "+grid_list) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Setting plotting options mpres@mpMaxLatF = Max_Lat mpres@mpMinLatF = Min_Lat mpres@mpMinLonF = Min_Lon mpres@mpMaxLonF = Max_Lon ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;; find xy model coordinates corresponding to lon/lat opts = True loc = wrf_user_ll_to_ij(a,(/Min_Lon,Max_Lon/),(/Min_Lat,Max_Lat/),opts) ;locmin = wrf_user_ll_to_ij(a, llmin(0),llmin(1), res) ;locmax = wrf_user_ll_to_ij(a, llmax(0),llmax(1), res) ; xmax = locmax(0) -1 ; edit for indices ; xmin = locmin(0) -1 ; ymax = locmax(1) -1 ; ymin = locmin(1) -1 ; nlat=ymax-ymin+1 ; correct for 0's, find size of arrays ; nlon=xmax-xmin+1 ; print(nlat) ; print(nlon) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (species .eq. "o3") var = wrf_user_getvar(a,"o3",it) var = var*1000 var@units = "ppb" end if if (isStrSubset(species,"tot")) ; this extracts all aerosol species below in the 8 bins var = get_aerosol_full(a, species, it) ; including water/liquid end if ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; set up model data to use ; z_sub = new((/nalt,nlat,nlon/),float) ; declare z for subregion ; z_sub = z(it,:,ymin:ymax,xmin:xmax) ; extract z values for subregion ;; mps_sub = new((/nalt,nlat,nlon/),float) ; declare pressure for subregion ; mps_sub = mps(it,:,ymin:ymax,xmin:xmax) ; extract pressure values for subregion ; var_sub = new((/nalt,nlat,nlon/),float) ; declare var for subregion & time ; var_sub = var(it,:,:,:) ; extract var values for current time ; u_sub = new((/nalt,nlat,nlon/),float) ; declare wind for subregion & time ; u_sub = u(it,:,:,:) ; extract wind values for current time ; v_sub = new((/nalt,nlat,nlon/),float) ; declare wind for subregion & time ; v_sub = v(it,:,:,:) ; extract wind values for current time ;;; set plot height ;;; ; height = plotheight(it) ; pressure = plotpressure(it) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Set resources for plots do level = 0, 1 ; loop over levels print("Working on level: " +level) opts = res opts@cnLineColor = "Blue" opts@cnLinesOn = False opts@cnFillOn = True opts@cnFillMode = "AreaFill" opts@gsnSpreadColors = True ; set basic plotting parameters: display_level = level + 1 ; eta level opts@PlotLevelID = "Eta Level " + display_level ;;;;;;;;;;;;;;;; PLOTTING PLOTS ;;;;;;;;;;;;;;;;;;;;; if (species .eq. "o3") opts@cnLevelSelectionMode = "ExplicitLevels" opts@FieldTitle = "o3_concentration at date: " + times(it) contour_var = wrf_contour(a,wks,var(level,:,:),opts) plot = wrf_map_overlays(a,wks,(/contour_var/),pltres,mpres) print("ploting o3") delete(contour_var) end if print("__________________________________________") end do ; end of level loop end do ;time loop print("Successful completion of program") end do ; file loop end do ; species end ; END OF PROGRAM