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/wrf/WRFUserARW.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" begin WRITE_NETCDF = True ; Whether to write output to NetCDF file ; TEXT FILE CONVERTING ROTATED AND SPHERICAL LATS AND LONS filename = "hirham_rot_geo_koor.dat" coords = asciiread(filename,(/100,110,4/),"float") lat2d = coords(:,:,2) lon2d = coords(:,:,3) printMinMax(lon2d, True) ;---MODEL file srcFileName = "atm.1961-02.daily.nc" sfile = addfile(srcFileName,"r") Opt = True Opt@SrcTitle = "HIRHAM" ; optional Opt@WgtFileName = "HIRHAM_to_Rect.nc" Opt@ForceOverwrite = True SLP = sfile->mslp U = sfile->u10 V = sfile->v10 lon2d = where(lon2d.lt.0,360+lon2d,lon2d) SLP@lat2d = lat2d ; This information will be used by SLP@lon2d = lon2d ; ESMF_regrid for the source grid U@lat2d = lat2d U@lon2d = lon2d V@lat2d = lat2d V@lon2d = lon2d dims = dimsizes(lat2d) nlat = dims(0) nlon = dims(1) Opt@SrcFileName = "WRF_SCRIP.nc" ; Name of source and Opt@DstFileName = "Rectilinear.nc" ; destination files ;---Create the destination lat/lon grid printVarSummary(lat2d) lat = fspan( 53.8, 90,nlat) lon = fspan( 0,360,nlon) lat@units = "degrees_north" lon@units = "degrees_east" lat!0 = "lat" lon!0 = "lon" lat&lat = lat lon&lon = lon Opt@DstGridType = "rectilinear" Opt@DstGridLat = lat Opt@DstGridLon = lon Opt@InterpMethod = "bilinear" Opt@SrcRegional = True Opt@DstRegional = True Opt@Debug = True SLP_regrid = ESMF_regrid(SLP,Opt) ; Do the regridding for TMP ; ; The source and destination grid description files and ; weight file will be the same for the next call to ; ESMF_grid, so no need to regenerate them. ; Opt@SkipSrcGrid = True Opt@SkipDstGrid = True Opt@SkipWgtGen = True U_regrid = ESMF_regrid(U,Opt) ; Do the regridding for V_regrid = ESMF_regrid(V,Opt) ;---Reset 0 values to missing values. SLP_regrid@_FillValue = default_fillvalue(typeof(SLP_regrid)) U_regrid@_FillValue = default_fillvalue(typeof(U_regrid)) V_regrid@_FillValue = default_fillvalue(typeof(V_regrid)) SLP_regrid = where(SLP_regrid.eq.0.0,SLP_regrid@_FillValue,\ SLP_regrid) U_regrid = where(U_regrid.eq.0.0,U_regrid@_FillValue,\ U_regrid) V_regrid = where(V_regrid.eq.0.0,V_regrid@_FillValue,\ V_regrid) printVarSummary(SLP_regrid) printVarSummary(U_regrid) printVarSummary(V_regrid) ;---------------------------------------------------------------------- ; Plotting section ; ; This section creates filled contour plots of both the original ; data and the regridded data, and panels them. ;---------------------------------------------------------------------- wks_slp = gsn_open_wks("png","interpolate_slp") wks_u = gsn_open_wks("png","interpolate_u") wks_v = gsn_open_wks("png","interpolate_v") res = True res@gsnDraw = False res@gsnFrame = False res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@lbLabelBarOn = False res@cnLevelSelectionMode = "ManualLevels" res@gsnPolar = "NH" res@mpMinLatF = min(lat) res@gsnAddCyclic = False nrec = 0 dims_orig = tostring(dimsizes(SLP(nrec,:,:))) mnmxint_slp = nice_mnmxintvl( min(SLP), max(SLP), 18, False) mnmxint_u = nice_mnmxintvl( min(U), max(U), 18, False) mnmxint_v = nice_mnmxintvl( min(V), max(V), 18, False) ;---SLP res@cnMinLevelValF = mnmxint_slp(0) res@cnMaxLevelValF = mnmxint_slp(1) res@cnLevelSpacingF = mnmxint_slp(2)/2. ; Create more levels res@tiMainFontHeightF = 0.015 ;---SLP regridded res@tiMainString = "rectilinear grid (" + Opt@InterpMethod + " interpolation)" slp_regrid = gsn_csm_contour_map(wks_slp,SLP_regrid(nrec,:,:),res) ;---SLP original res@tiMainString = "SLP on original curvilinear grid (" + \ str_join(dims_orig," x ") + ")" slp_orig = gsn_csm_contour_map(wks_slp,SLP(nrec,:,:),res) ;---U res@cnMinLevelValF = mnmxint_u(0) res@cnMaxLevelValF = mnmxint_u(1) res@cnLevelSpacingF = mnmxint_u(2)/2. ; Create more levels ;---U regridded res@tiMainString = "rectilinear grid (" + Opt@InterpMethod + " interpolation)" u_regrid = gsn_csm_contour_map(wks_u, U_regrid(nrec,:,:),res) ;---U original res@tiMainString = "U on original curvilinear grid (" + \ str_join(dims_orig," x ") + ")" u_orig = gsn_csm_contour_map(wks_u,U(nrec,:,:),res) ;---V plot res@cnMinLevelValF = mnmxint_v(0) res@cnMaxLevelValF = mnmxint_v(1) res@cnLevelSpacingF = mnmxint_v(2)/2. ; Create more levels res@tiMainString = "V on original curvilinear grid (" + \ str_join(dims_orig," x ") + ")" v_orig = gsn_csm_contour_map(wks_v,V(nrec,:,:),res) res@tiMainString = "rectilinear grid (" + Opt@InterpMethod + " interpolation)" v_regrid = gsn_csm_contour_map(wks_v, V_regrid(nrec,:,:),res) ;---Compare the plots in a panel pres = True pres@gsnMaximize = True pres@gsnPanelLabelBar = True gsn_panel(wks_slp,(/slp_orig,slp_regrid/),(/1,2/),pres) gsn_panel(wks_u,(/u_orig,u_regrid/),(/1,2/),pres) gsn_panel(wks_v,(/v_orig,v_regrid/),(/1,2/),pres) ;---Trim system("convert -trim interpolate_slp.png interpolate.png") system("convert -trim interpolate_u.png interpolate_u.png") system("convert -trim interpolate_v.png interpolate_v.png") if(WRITE_NETCDF) then dims_regrid = dimsizes(SLP_regrid) nrec = dims_regrid(0) nlat = dims_regrid(1) nlon = dims_regrid(2) fout_name = "example.nc" ; Output file system("rm -f " + fout_name) ; remove if exists fout = addfile (fout_name, "c") ; open output file ;=================================================================== ; explicitly declare file definition mode. Improve efficiency. ;=================================================================== setfileoption(fout,"DefineMode",True) ;=================================================================== ; create global attributes of the file ;=================================================================== fAtt = True ; assign file attributes fAtt@title = "NCL Efficient Approach to netCDF Creation" fAtt@source_file = srcFileName fAtt@description = "This file contains data regridded from the HIRHAM rotated lat-lon grid to a rectilinear grid" fAtt@creation_date = systemfunc ("date") fileattdef( fout, fAtt ) ; copy file attributes ;=================================================================== ; predefine the coordinate variables and their dimensionality ; Note: to get an UNLIMITED record dimension, we set the dimensionality ; to -1 (or the actual size) and set the dimension name to True. ;=================================================================== dimNames = getvardims(SLP_regrid) dimSizes = (/ -1, nlat, nlon/) dimUnlim = (/ True , False, False/) filedimdef(fout,dimNames,dimSizes,dimUnlim) ;=================================================================== ; predefine the the dimensionality of the variables to be written out ;=================================================================== ;=================================================================== filevardef(fout, "lat" ,typeof(lat),getvardims(lat)) filevardef(fout, "lon" ,typeof(lon),getvardims(lon)) filevardef(fout, "SLP" ,typeof(SLP_regrid),getvardims(SLP_regrid)) filevardef(fout, "U" ,typeof(U_regrid) ,getvardims(U)) filevardef(fout, "V" ,typeof(V_regrid) ,getvardims(V)) ;=================================================================== ; Copy attributes associated with each variable to the file ; All attributes associated with each variable will be copied. ;==================================================================== filevarattdef(fout,"SLP",SLP_regrid) ; copy SLP_regrid attributes filevarattdef(fout,"U", U_regrid) ; copy U_regrid attributes filevarattdef(fout,"V", V_regrid) ; copy V_regrid attributes filevarattdef(fout,"lat" ,lat) ; copy lat attributes filevarattdef(fout,"lon" ,lon) ; copy lon attributes ;=================================================================== ; explicitly exit file definition mode. **NOT REQUIRED** ;=================================================================== setfileoption(fout,"DefineMode",False) ;=================================================================== ; output only the data values since the dimensionality and such have ; been predefined. The "(/", "/)" syntax tells NCL to only output the ; data values to the predefined locations on the file. ;==================================================================== fout->lat = (/lat/) fout->lon = (/lon/) fout->SLP = (/SLP_regrid/) fout->U = (/U_regrid/) fout->V = (/V_regrid/) end if end