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" diri = "./" in = addfile(diri+fili,"r") ;names = getfilevarnames(in) ; Get all variable names on file ;print(names) time = in->initial_time0_hours ntim = dimsizes(time) date = cd_calendar(time, -3) ; yyyymmddhh lev = in->lv_ISBL1 ; 500, 700, 850, 1000 klev = dimsizes(lev) print(lev) ;************************************************ ; This high resolution grid is (800,1600) ;************************************************ nlat = 800 mlon = 2*nlat lat = latGau (nlat, "lat", "latitude", "degrees_north") lon = lonGlobeF (mlon, "lon", "longitude", "degrees_east") ;************************************************ ; loop over times & levels ; . use NCL synthesis function derive the variable ;************************************************ RH= new((/ntim,klev,nlat,mlon/), "float", 1e20) do nt=0,ntim-1 do kl=0,klev-1 RHC = in->R_GDS50_ISBL(nt,kl,:,:,:) ; complex ( 2, 800, 800) if (kt.eq.0 .and. kl.eq.0) then printVarSummary(RHC) end if RH(nt,kl,:,:) = shsgC(RHC, mlon ) ; synthesize and store end do ; kl end do ; nt printVarSummary(RH) ; raw array printMinMax(RH,0) ;************************************************ ; add meta data ;************************************************ RH!0 = "time" RH!1 = "lev" RH!2 = "lat" RH!3 = "lon" RH&time = time RH&lev = lev RH&lat = lat RH&lon = lon RH@long_name = RHC@long_name RH@units = RHC@units printVarSummary(RH) ; meta array printMinMax(RH,0) ;************************************************ ; compute a vertical average ;************************************************ ptop = lev(0) ; or set to 450 pbot = lev(klev-1) ; or set to 1010 dp = dpres_plevel(lev, pbot, ptop, 0) ; weights dp@longName = "layer thickness" dp@units = lev@units print(dp) ; layer thickness RHdp = dim_sum_wgt_n_Wrap(RH,dp,0,1) ; RH*dp printVarSummary(RHdp) RHavg = RHdp/sum(dp) ; vertical average copy_VarMeta(RHdp, RHavg) RHavg@long_name = "Vertical Avg: RH" ; add new attribute printVarSummary(RHavg) printMinMax(RHavg,0) ;************************************************ ; open workstation ;************************************************ wks = gsn_open_wks("ps", "rh.800x1600") ; open graphic wkstation gsn_define_colormap(wks,"BlAqGrYeOrRe") ; choose colormap res = True ; plot mods desired res@gsnSpreadColors = True ; use full range of color map res@gsnMaximize = True ; make as large as possible res@gsnPaperOrientation = "portrait" ; force portrait res@cnFillOn = True ; turn on color fill res@cnLinesOn = False ; turn off contour lines res@lbLabelAutoStride = True res@gsnCenterString = "" do nt=0,ntim-1 res@gsnCenterString = date(nt) plot = gsn_csm_contour_map_ce(wks,Tavg(nt,:,:),res) do kl=0,klev-1 res@gsnCenterString = lev(kl)+"hPa: "+date(nt) plot = gsn_csm_contour_map_ce(wks,RH(nt,kl,:,:),res) end do ; kl end do ; nt