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" load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" begin ;******************** READ WRF ****************************************** ;WRf files all_files = systemfunc ("ls /data14a/maurom/D_17/domain2/wrfout_d01_2017*") ;print(all_files) ;exit simulation = "D17_d02" longname = "Temperature_2m" compound = "T2" ; Remember to change the variables ncar_rea = "/data14a/maurom/D_17/meteo_MERRA_data/day_meteo_ncep/air.2m.gauss.2017.nc" title = compound title1 = compound + "_" + simulation diro = title1 system("mkdir -p " + diro) ;************************************************************************ ; The WRF ARW input file. a = addfiles (all_files, "r") ; Check Variables ;print(a) ;print(b1) ;exit ;************************************************************************ ; Variables WRF wr = wrf_user_getvar(a,compound,-1) wr = wr - 273.15 ; printVarSummary(wr) ;************************************************************************ b1 = addfile(ncar_rea,"r") lat = b1->lat ;;---Change (likely) lon = b1->lon ;;---Change (likely) ;print(b1) ;printVarSummary(src_lat) ;printVarSummary(src_lon) ;exit nc = b1->air(:,:,:) nc = nc - 273.15 nc = lonFlip(nc) ; change longitude (0;360) -> (-180;180) ; printVarSummary(nc) ;exit b = addfile ("/data14a/maurom/D_17/domain2/wrfout_d01_2017-01-01_00:00:00","r") dst_lat = b->XLAT(0,:,:) ;;---Change (maybe) dst_lon = b->XLONG(0,:,:) ;;---Change (maybe) printVarSummary(dst_lat) printVarSummary(dst_lon) ;exit ;************************************************************************ ;Times Variable 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 ;************************************************************************ ; Average wr&Time = times wr_avg = calculate_monthly_values(wr, "avg", 0, False) wr_avg@long_name= "WRFChem" wr_avg@units = "[degC]" delete(wr) nc_avg = calculate_monthly_values(nc, "avg", 0, False) nc_avg@long_name = "NCAR/NCEP" nc_avg@units = "[degC]" delete(nc) ; Check Variables ; printVarSummary(wr_avg) ; printVarSummary(nc_avg) ; 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 = wrf_map_resources(b,res) ; Add necessary resources for WRF map res@gsnDraw = False ; Don't draw plots res@gsnFrame = False ; Don't advance frame 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 = 35 res@cnLevelSpacingF = 5 ; 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 ;************************************************************************ ;---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 wr_avg&Time = cd_convert(wr_avg&Time,nc_avg&time@units) ;---Find the intersection between the two time arrays time_intersect = venn2_intersection(wr_avg&Time,nc_avg&time) copy_VarAtts(nc_avg&time,time_intersect) printVarSummary(time_intersect) print(time_intersect) ; 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 = wr_avg({time_intersect},:,:) nc_avg_subset = nc_avg({time_intersect},:,:) nc_regrid = rgrid2rcm_Wrap(lat, lon, nc_avg_subset, dst_lat, dst_lon, 0) ;copy_VarCoords(wr_avg_subset, nc_regrid) copy_VarMeta(wr_avg_subset, nc_regrid) printVarSummary(wr_avg_subset) printVarSummary(nc_avg_subset) printVarSummary(nc_regrid) ;exit 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_wrf = gsn_csm_contour_map(wks,wr_avg_subset(it,:,:),res_wrf) plot_nc = gsn_csm_contour_map(wks,nc_regrid(it,:,:),res_nc) ;---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 = "Monthly Average Temperature at 2m "+yyyymmddhh(it) gsn_panel(wks,(/plot_nc,plot_wrf/),(/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