;******************************************************************** ; rcm_7.ncl ; ; Concepts illustrated: ; - Read a specified variable ; - Extract date information via 'cd_calendar' ; - Explore the variable's statistical distribution ; - Determine the index corresponding to a user specified 'time' ; - Plotting RCM 'dust' data ; - Overlay wind onto plot ; - Extract required map parameters from the input file's attributes ; - Drawing color-filled contours over a lambert conformal map ;******************************************************************** ; float emflx(time, iy, jx) ; ; emflx:long_name = "Surface emission flux" ; ; emflx:standard_name = "surface_emission_flux" ; ; emflx:units = "mg m-2 day-1" ; ; =====> emflx:coordinates = "xlat xlon" ; ; =====> emflx:grid_mapping = "crs" ; ; emflx:cell_methods = "time: mean" ; ; ; The emflx:coordinates attribute indicates the appropriate lat/lon arrays ;******************************************************************** ; global attributes: ; :title = "ICTP Regional Climatic model V4" ; ; :institution = "ICTP" ; ; :source = "RegCM Model output file" ; ; :Conventions = "CF-1.4" ; ; :references = "http://gforge.ictp.it/gf/project/regcm" ; ; :model_revision = "tag-4.6.0" ; ; :history = "2019-08-12 13:49:21 : Created by RegCM RegCM Model program" ; ; :experiment = "Zabol-All-2018" ; ; =====> :projection = "LAMCON" ; ; :grid_size_in_meters = 20000. ; ; =====> :latitude_of_projection_origin = 31.025 ; ; =====> :longitude_of_projection_origin = 61.501 ; ; =====> :standard_parallel = 20., 44. ; ; [SNIP] ;******************************************************************** ; MAIN ;******************************************************************** d = "./" ; directory where files reside a = addfile(d+"Haboob-May-28-first_DUST04.2018052400.nc","r") b = addfile(d+"Haboob-May-28-first_ATM.2018052400.nc","r") ; Read in U and V at 1000 mb [subscript 0] ; x = a->emflx u = b->ua(:,10,:,:) v = b->va(:,10,:,:) printVarSummary(x) ; (time, iy, jx) printMinMax(x,0) print("-----") xlat = a->xlat ; (iy, jx) ==> Cartesian xlon = a->xlon dimx = dimsizes(x) ntim = dimx(0) nlat = dimx(1) mlon = dimx(2) ;---Extract dates: yyyymmddhh ymdh = cd_calendar(x&time, -3) print(ymdh) ymdhuv =cd_calendar(u&time, -3) print(ymdhuv) ; Explore data: look at distribution opt_sd = True opt_sd@PrintStat = True stat_x = stat_dispersion(x, opt_sd ) print("-----") ;---PLOT: Use "Cartesian" coordinates 'iy', 'ix' ; Use contour levels based upon data distribution. ; This may require iteration to find what is best. YMDH = 2018052812 ; desired date nt = ind(YMDH.eq.ymdh) ; time index corresponding to the requested data ut = ind(YMDH.eq.ymdhuv) pltTyp = "png" ; png, x11, pdf, ... pltDir = "./" ;;pltFil = "rcm_"+YMDH pltFil = "rcm" pltPth = pltDir+pltFil ; Delete cartesian coordinate arrays associated with the variables ; For this application, NCL does not use this information delete( [/ x&jx, x&iy, u&jx, u&iy, v&jx, v&iy /] ) wks = gsn_open_wks(pltTyp,pltPth) res = True ; plot mods desired res@gsnMaximize = True ; maximize plot size res@gsnAddCyclic = False ; data are regional res@cnFillOn = False ; turn on color res@cnFillOn = True ; turn on color fill res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off contour line labels res@cnFillPalette = "WhiteBlueGreenYellowRed" res@cnFillMode = "RasterFill" ; Raster Mode ;;res@lbOrientation = "Vertical" res@lbBoxLinesOn = False ; Turn off labelbar box lines res@mpFillOn = False res@mpLimitMode = "Corners" ; choose region of map res@mpLeftCornerLatF = xlat(0,0) res@mpLeftCornerLonF = xlon(0,0) res@mpRightCornerLatF = xlat(nlat-1,mlon-1) res@mpRightCornerLonF = xlon(nlat-1,mlon-1) res@mpDataBaseVersion = "MediumRes" ; better map outlines res@pmTickMarkDisplayMode= "Always" ; turn on tickmarks ; Map information parameters are file attributes [also associated with 'crs' variable res@mpProjection = "LambertConformal" res@mpLambertParallel1F = a@standard_parallel(0) res@mpLambertParallel2F = a@standard_parallel(1) res@mpLambertMeridianF = a@longitude_of_projection_origin res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/ 100, 250, 500, 750, 1000, 2000\ , 3000, 4000, 5000, 6000, 7000, 8000\ , 9000,10000,11000,12000,13000,14000\ , 15000,17250,20000,25000,30000,40000\ , 50000,60000,75000,100000 /) res@sfXArray = xlon ; equivalently x@lon2d = xlon res@sfYArray = xlat ; equivalently x@lat2d = xlat res@tiMainString = "ITCP: RCM: "+YMDH res@gsnDraw = False ; turn off draw and frame res@gsnFrame = False ; b/c this is an overlay plot plotc = gsn_csm_contour_map(wks,x(nt,:,:),res) resv = True resv@gsnDraw = False ; turn off draw and frame resv@gsnFrame = False ; b/c this is an overlay plot resv@vcRefMagnitudeF = 5.0 ; define vector ref mag resv@vcRefLengthF = 0.045 ; define length of vec ref resv@vcRefAnnoOrthogonalPosF = -0.5 ; move ref vector resv@vcRefAnnoArrowLineColor = "black" ; change ref vector color resv@vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref resv@vcGlyphStyle = "CurlyVector" ; turn on curly vectors ;resv@vcLineArrowColor = "white" ; change vector color resv@vcLineArrowColor = "black" ; change vector color ;resv@vcLineArrowThicknessF = 2.0 ; change vector thickness ;resv@vcVectorDrawOrder = "PostDraw" plotv=gsn_vector(wks,u(ut,:,:),v(ut,:,:),resv) overlay(plotc,plotv) draw(plotc) frame(wks)