;************************************************************************************** ; Script produced by Justin Glisan at Iowa State University for the purpose of cutting ; out the full SN extreme analysis box (Over North America). After this is completed, ; masking commands are used to mask the the ocean and sea ice from Pan-Arctic WRF data. ; The data that is left is strictly land-based precipitation. ;*************************************************************************************** ;*********************************************** 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/wrf/WRF_contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;load "/usr/local/ncl/lib/ncarg/nclscripts/csm/shea_util.ncl" ;load "/usr/local/ncl/lib/ncarg/nclscripts/csm/contributed.ncl" ;************************************************ begin ;************************************************ ; create pointer to file and read in data ;************************************************ a = addfile("/archive/u1/uaf/glisan/CORDEX_OUTPUT/CA_Production/Precip_Daily/wrfout_CA_TP_F_JAS_16.nc","r") m = addfile("/archive/u1/uaf/glisan/CORDEX_OUTPUT/A_Ensemble_Pingo/1989A/wrfpost_LANDMASK_198901.nc","r") tt = addfile("/archive/u1/uaf/glisan/CORDEX_OUTPUT/SEASON/A/wrfpost_A_JAS_16.nc","r") wks = gsn_open_wks("pdf","CORDEX_ARCTIC_Mask_test") ;************************************************ ; Set some universal time and lat/lon. coords ;************************************************ t = 5823; Ensemble member= t - 4 TimeChara = tt->time(0:5823) DimTimeChara = dimsizes(TimeChara) ; print("TimeChara dimsizes:" + DimTimeChara) nTime = TimeChara nTimes = dimsizes(nTime) - 1 ; print("nTime:" + nTime + "and nt: " + nTimes) lat = m->lat lon = m->lon lat2d = lat lon2d = lon res = True ;************************************************ ; Read in precip. data - three "variables since ; the largest sub-region needs three boxes to ; get produce the proper cut-out. ;************************************************ wrf_mean = a->PAW_TP_Mean(0:t,:,:) wrf_mean!0 = "time" wrf_mean!1 = "lat" wrf_mean!2 = "lon" wrf_mean1 = a->PAW_TP_Mean(0:t,:,:) wrf_mean1!0 = "time" wrf_mean1!1 = "lat" wrf_mean1!2 = "lon" wrf_mean2 = a->PAW_TP_Mean(0:t,:,:) wrf_mean2!0 = "time" wrf_mean2!1 = "lat" wrf_mean2!2 = "lon" ;************************************************ ; Define the new arrays so daily averages ; can be computed over the 16 year period. ;************************************************ Daily_2mT_PAW = new ( (/1456,135,125/), "float") Daily_2mT_PAW1 = new ( (/1456,135,125/), "float") Daily_2mT_PAW1 = new ( (/1456,135,125/), "float") ;************************************************ ; Averages are taken here...data every 6-hrs, ; so four t.s. averages for daily ;************************************************ Daily_2mT_PAW = wrf_mean(::4,:,:) ntJump = 4 ; 6-hr step ntStrt = 0 ; start at 0 time ntLast = ntJump-1 ; end every 4 steps or 1 day do n=0,nTimes,4 Daily_2mT_PAW(n/4,:,:) = dim_avg_Wrap(wrf_mean(lat|:,lon|:,time|ntStrt:ntLast)) ntStrt = ntStrt+ntJump ntLast = ntLast+ntJump end do delete(ntStrt) delete(ntLast) delete(ntJump) delete(n) Daily_2mT_PAW1 = wrf_mean1(::4,:,:) ntJump = 4 ; 6-hr step ntStrt = 0 ; start at 0 time ntLast = ntJump-1 ; end every 4 steps or 1 day do n=0,nTimes,4 Daily_2mT_PAW1(n/4,:,:) = dim_avg_Wrap(wrf_mean1(lat|:,lon|:,time|ntStrt:ntLast)) ntStrt = ntStrt+ntJump ntLast = ntLast+ntJump end do delete(ntStrt) delete(ntLast) delete(ntJump) delete(n) Daily_2mT_PAW2 = wrf_mean2(::4,:,:) ntJump = 4 ; 6-hr step ntStrt = 0 ; start at 0 time ntLast = ntJump-1 ; end every 4 steps or 1 day do n=0,nTimes,4 Daily_2mT_PAW2(n/4,:,:) = dim_avg_Wrap(wrf_mean2(lat|:,lon|:,time|ntStrt:ntLast)) ntStrt = ntStrt+ntJump ntLast = ntLast+ntJump end do delete(ntStrt) delete(ntLast) delete(ntJump) delete(n) ;********************************************** ; The next six lines are added to streamline ; an existing script (quick and dirty fix) ;********************************************** delete(wrf_mean) delete(wrf_mean1) delete(wrf_mean2) wrf_mean = Daily_2mT_PAW wrf_mean1 = Daily_2mT_PAW1 wrf_mean2 = Daily_2mT_PAW2 ;************************************************ ; Read seaice, land, and ocean flags to creates ; the approp. masks ;************************************************ seaice = m->SeaIce seaice!0 = "time" seaice!1 = "lat" seaice!2 = "lon" seaice1 = m->SeaIce seaice1!0 = "time" seaice1!1 = "lat" seaice1!2 = "lon" seaice2 = m->SeaIce seaice2!0 = "time" seaice2!1 = "lat" seaice2!2 = "lon" land = m->LandMask land!0 = "time" land!1 = "lat" land!2 = "lon" land1 = m->LandMask land1!0 = "time" land1!1 = "lat" land1!2 = "lon" land2 = m->LandMask land2!0 = "time" land2!1 = "lat" land2!2 = "lon" ;************************************************ ; Time average "flag" data for masks ;************************************************ land_mean = dim_avg( land(lat|:, lon|:, time|:) ) land_mean1 = dim_avg( land1(lat|:, lon|:, time|:) ) land_mean2 = dim_avg( land2(lat|:, lon|:, time|:) ) seaice_mean = dim_avg( seaice(lat|:, lon|:, time|:) ) seaice_mean1 = dim_avg( seaice1(lat|:, lon|:, time|:) ) seaice_mean2 = dim_avg( seaice2(lat|:, lon|:, time|:) ) ;************************************************ ; Set new attributes for averages ;************************************************ wrf_mean@lat2d = lat ; 2D coordinate arrays wrf_mean@lon2d = lon wrf_mean1@lat2d = lat ; 2D coordinate arrays wrf_mean1@lon2d = lon wrf_mean2@lat2d = lat ; 2D coordinate arrays wrf_mean2@lon2d = lon land_mean@lat2d = lat ; 2D coordinate arrays land_mean@lon2d = lon land_mean1@lat2d = lat ; 2D coordinate arrays land_mean1@lon2d = lon land_mean2@lat2d = lat ; 2D coordinate arrays land_mean2@lon2d = lon seaice_mean@lat2d = lat ; 2D coordinate arrays seaice_mean@lon2d = lon seaice_mean1@lat2d = lat ; 2D coordinate arrays seaice_mean1@lon2d = lon seaice_mean2@lat2d = lat ; 2D coordinate arrays seaice_mean2@lon2d = lon ;************************************************ ; Here, the NA domain is selected. Within this ; full domain, the four sub-boxes are found. ; For each box, input selected lat/lon for ; the necessary box(es). ;************************************************ latMin = 70 latMax = 85 latMin1 = 68 latMax1 = 85 latMin2 = 50 latMax2 = 68 lonMin = -126 lonMax = -119 lonMin1 = -119 lonMax1 = -70 lonMin2 = -107 lonMax2 = -53 ;************************************************ ; Next, we slice the domain out of each field ;************************************************ wrf_mean = mask(wrf_mean,(wrf_mean@lat2d.ge.latMin .and. wrf_mean@lat2d.le.latMax .and. \ wrf_mean@lon2d.ge.lonMin .and. wrf_mean@lon2d.le.lonMax), True) wrf_mean1 = mask(wrf_mean1,(wrf_mean1@lat2d.ge.latMin1 .and. wrf_mean1@lat2d.le.latMax1 .and. \ wrf_mean1@lon2d.ge.lonMin1 .and. wrf_mean1@lon2d.le.lonMax1), True) wrf_mean2 = mask(wrf_mean2,(wrf_mean2@lat2d.ge.latMin2 .and. wrf_mean2@lat2d.le.latMax2 .and. \ wrf_mean2@lon2d.ge.lonMin2 .and. wrf_mean2@lon2d.le.lonMax2), True) seaice_mean = mask(seaice_mean,(seaice_mean@lat2d.ge.latMin .and. seaice_mean@lat2d.le.latMax .and. \ seaice_mean@lon2d.ge.lonMin .and. seaice_mean@lon2d.le.lonMax), True) seaice_mean1 = mask(seaice_mean1,(seaice_mean1@lat2d.ge.latMin1 .and. seaice_mean1@lat2d.le.latMax1 .and. \ seaice_mean1@lon2d.ge.lonMin1 .and. seaice_mean1@lon2d.le.lonMax1), True) seaice_mean2 = mask(seaice_mean2,(seaice_mean2@lat2d.ge.latMin2 .and. seaice_mean2@lat2d.le.latMax2 .and. \ seaice_mean2@lon2d.ge.lonMin2 .and. seaice_mean2@lon2d.le.lonMax2), True) land_mean = mask(land_mean,(land_mean@lat2d.ge.latMin .and. land_mean@lat2d.le.latMax .and. \ land_mean@lon2d.ge.lonMin .and. land_mean@lon2d.le.lonMax), True) land_mean1 = mask(land_mean1,(land_mean1@lat2d.ge.latMin1 .and. land_mean1@lat2d.le.latMax1 .and. \ land_mean1@lon2d.ge.lonMin1 .and. land_mean1@lon2d.le.lonMax1), True) land_mean2 = mask(land_mean2,(land_mean2@lat2d.ge.latMin2 .and. land_mean2@lat2d.le.latMax2 .and. \ land_mean2@lon2d.ge.lonMin2 .and. land_mean2@lon2d.le.lonMax2), True) ;************************************************ ; use mask function to mask out seaice (seaice=0 ; from ocean data and then ocean data from land ; ocean = 1 ;************************************************ seaice_only = mask(wrf_mean,seaice_mean,0) seaice_only@lat2d = lat ; 2D coordinate arrays seaice_only@lon2d = lon seaice_only1 = mask(wrf_mean1,seaice_mean1,0) seaice_only1@lat2d = lat ; 2D coordinate arrays seaice_only1@lon2d = lon seaice_only2 = mask(wrf_mean2,seaice_mean2,0) seaice_only2@lat2d = lat ; 2D coordinate arrays seaice_only2@lon2d = lon land_mask = mask(seaice_only,land_mean,1) land_mask@lat2d = lat ; 2D coordinate arrays land_mask@lon2d = lon land_mask!0 = "time" land_mask!1 = "lat" land_mask!2 = "lon" land_mask1 = mask(seaice_only1,land_mean1,1) land_mask1@lat2d = lat ; 2D coordinate arrays land_mask1@lon2d = lon land_mask1!0 = "time" land_mask1!1 = "lat" land_mask1!2 = "lon" land_mask2 = mask(seaice_only2,land_mean2,1) land_mask2@lat2d = lat ; 2D coordinate arrays land_mask2@lon2d = lon land_mask2!0 = "time" land_mask2!1 = "lat" land_mask2!2 = "lon" ;**************************************************************** ; Use the x_crop command to "loop" around the selected ; "square" boundaries for each box. Since everything outside ; of the sub-regional boxes will be outputted as missing values, ; these commands will effectively parse the data by cropping the ; full CORDEX Arctic domain to the closest grid squares to the ; boxes. Some missing values still outputted, but a vast majority ; are not! ;***************************************************************** i = ind_resolve(ind(.not.ismissing(ndtooned(land_mask))),dimsizes(land_mask)) x_crop = land_mask( min(i(:,0)):max(i(:,0)), min(i(:,1)):max(i(:,1)), min(i(:,2)):max(i(:,2)) ) x_crop@lat2d = lat ; 2D coordinate arrays x_crop@lon2d = lon x_crop!0 = "time" x_crop!1 = "lat" x_crop!2 = "lon" ii = ind_resolve(ind(.not.ismissing(ndtooned(land_mask1))),dimsizes(land_mask1)) x_crop1 = land_mask1( min(ii(:,0)):max(ii(:,0)), min(ii(:,1)):max(ii(:,1)), min(ii(:,2)):max(ii(:,2)) ) x_crop1@lat2d = lat ; 2D coordinate arrays x_crop1@lon2d = lon x_crop1!0 = "time" x_crop1!1 = "lat" x_crop1!2 = "lon" iii = ind_resolve(ind(.not.ismissing(ndtooned(land_mask2))),dimsizes(land_mask2)) x_crop2 = land_mask2( min(iii(:,0)):max(iii(:,0)), min(iii(:,1)):max(iii(:,1)), min(iii(:,2)):max(iii(:,2)) ) x_crop2@lat2d = lat ; 2D coordinate arrays x_crop2@lon2d = lon x_crop2!0 = "time" x_crop2!1 = "lat" x_crop2!2 = "lon" print(x_crop) print(x_crop1) print(x_crop2) ;******************************************************** ; Have not automated script, yet, so I print the box(es) ; using the simple print command and just push the output ; into text files. ;******************************************************** ;************************************************ ; Here plotting attributes are set ;************************************************ res = True res@tiMainFontHeightF = 0.012 res@tiMainPosition = "Left" res@pmTickMarkDisplayMode = "Never" res@pmTickMarkZone = 0.1 res@pmLabelBarOrthogonalPosF = 0.10 res@tmXTOn = False res@tmXBOn = False res@tmYROn = False res@tmXTOn = False res@gsnMaximize = True res@lbOrientation = "Vertical" res@lbTitlePosition = "Right" res@lbTitleDirection = "Across" res@lbTitleAngleF = 90 res@lbTitleFontHeightF = 0.015 res@lbTitleString = "mm/month" res@lbTitleOffsetF = -0.001 res@lbBoxMinorExtentF = 0.12 ;.20 res@lbLabelPosition = "Right" res@lbLabelAutoStride = True res@lbLabelDirection = "Across" res@lbLabelAngleF = 0 res@mpProjection = "Stereographic" res@mpProjection = "Stereographic" res@mpCenterLatF = 88.86 res@mpCenterLonF = 179.4 ; allignes the plot squarely res@mpGridSpacingF = 10 res@mpGridPolarLonSpacingF = 20 res@mpLimitMode = "NPC" res@mpLeftNPCF = 0.3525 res@mpBottomNPCF = 0.365 res@mpRightNPCF = 0.753 res@mpTopNPCF = 0.647 res@mpLimitMode = "Corners" res@mpLeftCornerLonF = -44.9 res@mpLeftCornerLatF = 51.9 res@mpRightCornerLatF = 50.1 res@mpRightCornerLonF = 138.5 res@mpFillOn = False res@mpGridAndLimbOn = True res@mpOutlineBoundarySets = "AllBoundaries" res@gsnMaximize = True ; res@gsnPolar = "NH" res@gsnFrame = True res@mpFillOn = False res@mpGridAndLimbOn = True res@mpOutlineBoundarySets = "AllBoundaries" res@gsnMaximize = True ; res@gsnPolar = "NH" res@gsnFrame = True res@mpFillPatternBackground = -1 res@mpGeophysicalLineColor = "Black" res@mpGeophysicalLineThicknessF = 1.0 res@mpGridAndLimbOn = True res@mpGridMaskMode = "MaskLand" res@mpGridLineDashPattern = 2 res@mpPerimOn = True res@mpPerimLineColor = "Black" res@mpPerimLineThicknessF = 1.5 res@gsnDraw = False res@gsnFrame = False gsn_define_colormap(wks,"WhViBlGrYeOrRe") ; gsn_define_colormap(wks,"BlWhRe") ;************************************************* ;**CONFIGURE BELOW******************************** ;************************************************* res@gsnSpreadColors = True res@gsnSpreadColorEnd = -1 res@cnFillOn = True ; res@cnLevelSelectionMode = "ExplicitLevels" ; res@cnLevels = (/0,100,200,300,400,500,600,700,800/) res@cnLinesOn = False plot = gsn_csm_contour_map(wks,x_crop(70,:,:),res) txt0 = create "MainPlotTitle" textItemClass wks "txString" : "2007-07 Triple SN PAW Mean Total Accumulated Precip." "txFontHeightF" : 0.015 end create anno = NhlAddAnnotation(plot,txt0) setvalues anno "amZone" : 3 "amSide" : "Top" "amJust" : "CenterLeft" "amParallelPosF" : 0.0 "amOrthogonalPosF": 0.0 "amResizeNotify" : False end setvalues draw(plot) frame(wks) ;************************************************ ; Let's clean up and set output (if needed) ;************************************************ ; delete(res) ; delete(txt0) ; delete(path) ; fileAtt = True ; fileAtt@creation_date = systemfunc("date") ; fileAtt@created_by = "Justin Glisan - glisanj@iastate.edu" ; fileAtt@comment0 = "Created with NCL script: Ocean_Seaice_Mask.ncl" ; fileattdef(fout, fileAtt) end