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/esmf/ESMF_regridding.ncl" ;; Needed for regridding begin Opt = True Opt@ForceOverwrite = True Opt@DstGridType = "1.0x1.0" Opt@DstLLCorner = (/ -90, 0 /) ;; lower-left corner lat/lon Opt@DstURCorner = (/ 90, 359 /) ;; upper-right corner lat/lon latRange = (/0,90/) lonRange = (/0,359/) ;;; increments at 1.0 degree intervals domain_lat = fspan(latRange(0),latRange(1),91) domain_lon = fspan(lonRange(0),lonRange(1),360) ;; sort of a dummy array to get the 2-d grid count_total = new((/dimsizes(domain_lat),dimsizes(domain_lon)/),"integer") count_total!0 = "lat" count_total!1 = "lon" count_total&lat = domain_lat count_total&lon = domain_lon count_total = 0 ;; lats and lons of domain, dimensioned same as count_total lat_tot = conform(count_total,count_total&lat,0) lon_tot = conform(count_total,count_total&lon,1) rad = 4.0*atan(1.0)/180.0 clat = cos(domain_lat*rad) fileName_skill = "/keyserlab_rit/kbiernat/arctic_forecast_skill/GEFS_JJA/low_high_skill_aa_spread_aa_rmse_v2_files/low_skill_periods_10per_aa_rmse_v2_JJA_v2.nc" inFile_skill = addfile( fileName_skill, "r" ) time_skill = inFile_skill->time(:) ntim_skill = dimsizes(time_skill) inFile_ACs = addfile("/keyserlab_rit/kbiernat/arctic_cyclone_tracks/arctic_cyclones_JJA_skill/arctic_cyclones_low_skill_periods_10per_aa_rmse_v2_JJA_v2.nc","r") time_old = inFile_ACs->time(:,:) lat_old = inFile_ACs->lat(:,:) lon_old = inFile_ACs->lon(:,:) pmin_old = inFile_ACs->pmin(:,:) pcont_old = inFile_ACs->pcont(:,:) idtrack_old = inFile_ACs->idtrack(:,:) lenTrack_old = inFile_ACs->lenTrack(:) ;; 2d array where values are (nTracks,nTimesTrack) nTracks = dimsizes(lenTrack_old) nTimesMax = max(lenTrack_old) largest_aa_pos_div_avg = new((/nTracks/), float) lat = new((/nTracks/), float) lon = new((/nTracks/), float) pmin = new((/nTracks/), float) pcont = new((/nTracks/), float) idtrack = new((/nTracks/), float) time = new((/nTracks/), double) lenTrack = new((/nTracks/), integer) time@units = "hours since 1800-01-01 00:00:00" levs_div = (/350,300,250/) do i=0, nTracks-1 maxVal = -9999999.0 print("On track "+i) do j=0, lenTrack_old(i)-1 do k=0, ntim_skill-1 if (time_old(i,j).eq.time_skill(k)) then if (lat_old(i,j).gt.70.0) then timeCurr = time_old(i,j) timeCurr@units = "hours since 1800-01-01 00:00:00" utc_time := ut_calendar(timeCurr,0) yyyy = floattointeger(utc_time(0,0)) mm = floattointeger(utc_time(0,1)) dd = floattointeger(utc_time(0,2)) hh = floattointeger(utc_time(0,3)) u_file = addfile( "/erai/"+yyyy+"/u."+yyyy+".nc" , "r" ) v_file = addfile( "/erai/"+yyyy+"/v."+yyyy+".nc", "r" ) u = u_file->u({timeCurr},{levs_div},:,:) v = v_file->v({timeCurr},{levs_div},:,:) div = uv2dvG_Wrap(u,v) ; calculate divergence from the u and v wind components div_avg_old = dim_avg_n_Wrap(div,0) div_avg_old = div_avg_old* 10^5 div_avg = ESMF_regrid_with_weights(div_avg_old,"weights_file.nc",Opt) div_avg_domain := div_avg({domain_lat},{domain_lon}) pos_div_avg = div_avg_domain pos_div_avg = where(div_avg_domain.gt.0,div_avg_domain,div_avg_domain@_FillValue) lat_AC = count_total&lat lon_AC = count_total&lon lat_AC = lat_old(i,j) lon_AC = lon_old(i,j) ;;----------------------------------------------------------------- ;;Calculate Great Circle Distances between points ;;----------------------------------------------------------------- r = 6370. ; radius of spherical Earth (km) d2r = 0.0174532925199433 ; degrees to radians conversion factor lat1 = lat_tot lon1 = lon_tot lat2 = conform(count_total,lat_AC,0) lon2 = conform(count_total,lon_AC,1) lat1 = lat1*d2r lon1 = lon1*d2r lat2 = lat2*d2r lon2 = lon2*d2r dlat = lat2-lat1 dlon = lon2-lon1 latTerm = sin(.5*dlat) latTerm = latTerm*latTerm lonTerm = sin(.5*dlon) lonTerm = lonTerm*lonTerm*cos(lat1)*cos(lat2) dAngle = sqrt(latTerm+lonTerm) dist = 2.*r*asin(dAngle) ;;----------------------------------------------------------------- rad_search = 1000.0 ; radius we search outward to (km) pos_div_avg_inCircle = where((dist.lt.rad_search),pos_div_avg(:,:),pos_div_avg@_FillValue) pos_div_avg_inCircle@_FillValue = pos_div_avg@_FillValue aa_pos_div_avg = wgt_areaave(pos_div_avg_inCircle, clat, 1.0, 0) if( ismissing(aa_pos_div_avg)) then delete(u_file) delete(v_file) delete(u) delete(v) delete(div) delete(div_avg_old) delete(div_avg) delete(div_avg_domain) delete(pos_div_avg) delete(pos_div_avg_inCircle) delete(aa_pos_div_avg) continue end if if (aa_pos_div_avg.gt.maxVal) then largest_aa_pos_div_avg(i) = aa_pos_div_avg time(i) = time_old(i,j) lat(i) = lat_old(i,j) lon(i) = lon_old(i,j) pmin(i) = pmin_old(i,j) pcont(i) = pcont_old(i,j) idtrack(i) = idtrack_old(i,j) lenTrack(i) = lenTrack_old(i) maxVal = aa_pos_div_avg end if delete(u_file) delete(v_file) delete(u) delete(v) delete(div) delete(div_avg_old) delete(div_avg) delete(div_avg_domain) delete(pos_div_avg) delete(pos_div_avg_inCircle) delete(aa_pos_div_avg) end if end if end do end do end do ;;################################################## ;;Output variable to netCDF file ;;################################################## out_fil = "/keyserlab_rit/kbiernat/v11_1_Crawford_cyclone_tracking/cyclonetracking/gefs_JJA/largest_div_350_to_250hPa/files/largest_aa_pos_div_avg_350_to_250hPa_ACs_during_low_skill_periods_inArctic_v2.nc" ;-Output filename var1 = "largest_aa_pos_div_avg" var2 = "lat" var3 = "lon" var4 = "time" var5 = "pmin" var6 = "pcont" var7 = "idtrack" var8 = "lenTrack" ;--------- Output a compressed netCDF ------ setfileoption("nc","Format","NetCDF4Classic") setfileoption("nc","CompressionLevel",1) system("/bin/rm -f "+out_fil) fout = addfile(out_fil ,"c") dimNames = (/"nTracks"/) dimSizes = (/ nTracks /) dimUnlim = (/ False /) filedimdef(fout,dimNames,dimSizes,dimUnlim) chunks = (/ nTracks /) filechunkdimdef(fout,dimNames,chunks,dimUnlim) filevardef(fout, var1 ,typeof(largest_aa_pos_div_avg), (/"nTracks"/)) filevardef(fout, var2 ,typeof(lat),(/"nTracks"/)) filevardef(fout, var3 ,typeof(lon),(/"nTracks"/)) filevardef(fout, var4 ,typeof(time),(/"nTracks"/)) filevardef(fout, var5 ,typeof(pmin),(/"nTracks"/)) filevardef(fout, var6 ,typeof(pcont),(/"nTracks"/)) filevardef(fout, var7 ,typeof(idtrack),(/"nTracks"/)) filevardef(fout, var8 ,typeof(lenTrack),(/"nTracks"/)) filevarattdef(fout,"time",time) ; copy time attributes fout->$var1$(:) = (/largest_aa_pos_div_avg/) fout->$var2$(:) = (/lat/) fout->$var3$(:) = (/lon/) fout->$var4$(:) = (/time/) fout->$var5$(:) = (/pmin/) fout->$var6$(:) = (/pcont/) fout->$var7$(:) = (/idtrack/) fout->$var8$(:) = (/lenTrack/) delete(dimNames) delete(dimSizes) delete(dimUnlim) delete(chunks) end