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/csm/shea_util.ncl" external XMISS "./xmiss.so" ; this script takes the monthly averaged CAM data ; and computes seasonal averages begin ;path = "/scratch/scratch96/j/jtrapp/gcm-data/030bES_T85/" ;path = "/autohome/u115/cvillanu/Past/1970-1999_z_p/" path = "./" ;****************************************************** ; get some initial information from one of the netcdf files ;****************************************************** decade = 1980 ;filstrng = "030bES_T85" main1= "/autohome/u115/cvillanu/NARR/" ;main1 = "/scratch/scratch96/j/jtrapp/gcm-data/"+filstrng+"/" ;main1 = "/usr/rmt_share/scratch77/cecille/" file1=decade+"Jun-narr.nc" ;Pointers f = addfile(main1+file1, "r") ;Read Variables ... ; mainly to extract some metadata, to get dimensions lon = f->lon lat = f->lat ;lev = f->lev(4:25) lev = f->lv_ISBL3 tt = f->avgU3D nlat = dimsizes(lat) nlon = dimsizes(lon) nlev = dimsizes(lev) nz = nlev-1 print ("number of pressure levels="+nlev) print ("number of lat="+nlat) print ("number of lon="+nlon) monstrng = (/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/) ;Set Path ;t_mam = new((/nlev,nlat,nlon/), float) t_jja = new((/nlev,nlat,nlon/), float) ;t_jja = new((/nlon,nlat,nlev/), float) ;q_mam = new((/nlev,nlat,nlon/), float) q_jja = new((/nlev,nlat,nlon/), float) ;u_mam = new((/nlev,nlat,nlon/), float) u_jja = new((/nlev,nlat,nlon/), float) ;v_mam = new((/nlev,nlat,nlon/), float) v_jja = new((/nlev,nlat,nlon/), float) ;p_mam = new((/nlev,nlat,nlon/), float) ;p_jja = new((/nlev,nlat,nlon/), float) ;ps_mam = new((/nlat,nlon/), float) ps_jja = new((/nlat,nlon/), float) ;z_mam = new((/nlev,nlat,nlon/), float) z_jja = new((/nlev,nlat,nlon/), float) ;cape_mam = new((/nlat,nlon/), float) cape_jja = new((/nlat,nlon/), float) ;sev_mam = new((/nlat,nlon/), float) ;sev_jja = new((/nlat,nlon/), float) ;******************************************************************* ; these two arrays will count the number of times, at each gridpoint ; where there are *valid* (i.e., nonmissing) values ;******************************************************************* miss_jja = new((/nlat,nlon/), float) ;miss_mam = new((/nlat,nlon/), float) ;t_mam = 0. t_jja = 0. ;q_mam = 0. q_jja = 0. ;u_mam = 0. u_jja = 0. ;v_mam = 0. v_jja = 0. ;p_mam = 0. ;p_jja = 0. ;ps_mam = 0. ps_jja = 0. ;z_mam = 0. z_jja = 0. ;cape_mam = 0. cape_jja = 0. ;sev_mam = 0. ;sev_jja = 0. ;miss_mam = 0. miss_jja = 0. ;mamcnt = 0. jjacnt = 0. ;************************************************ ; create plot ;************************************************ wks = gsn_open_wks("pdf","means_NARR") ; open a ps file wks = gsn_open_wks("ncgm","means_NARR") ; open a ps file gsn_define_colormap (wks,"gui_default") ; choose color map res = True ; plot mods desired res@gsnSpreadColors = True res@cnFillOn = True ; turn on color fill res@mpFillOn = False ; turn off continent gray res@mpOutlineBoundarySets = "GeophysicalAndUSStates" res@cnLinesOn = True res@cnLineLabelsOn = True ; res@gsnAddCyclic = False ; this applies to the windowed domain res@mpLimitMode = "LatLon" ; control map zoom with lat/lons res@mpMaxLatF = 55. ; actual values of zoom res@mpMinLatF = 25. res@mpMinLonF = 220. res@mpMaxLonF = 300. res@gsnMaximize = False ; res@cnLevelSelectionMode = "ManualLevels" ; res@cnLevelSpacingF = 5. ; res@cnMinLevelValF = 1. ; res@cnMaxLevelValF = 50. ; ; loop through files and extract the monnthly data ;ybeg = 1978 ;yend = 1998 ybeg = 1980 yend = 1981 do ny = ybeg,yend ;do nm = 2, 4 ;March-August do nm = 5, 7 ;March-August ;do nm = 5, 5 ;March-August file3 = ny+monstrng(nm)+"-narr.nc" ;Pointer f3 = addfile(path+file3, "r") ; check for missing data (essentially, a missing file) if(any(.not.ismissing(f3))) then print ("file3="+file3) ; read your variables here: T, Q, U, V, Z, P ; these are dimensions: nlev, nlat, nlon t = f3->avgT3D q = f3->avgQ3D u = f3->avgU3D v = f3->avgV3D ; p = f3->pavg ps = f3->avgps z = f3->avgz ;sev = f3->sev cape = f3->avgtcape ; print("the missing value for t is"+t@_FillValue) ; plot = gsn_csm_contour_map(wks,sev, res) ; create plot ; plot = gsn_csm_contour_map(wks,z(nz,:,:), res) ; create plot ; if (nm.ge.2.and.nm.le.4) then ; choose any of the variables to determine the nonmissing values ; XMISS::xmiss(t(0,:,:), miss_mam, nlon, nlat) ; now, go through and set missing values to zero ; t=where(ismissing(t),0.,t) ; q=where(ismissing(q),0.,q) ; u=where(ismissing(u),0.,u) ; v=where(ismissing(v),0.,v) ; p=where(ismissing(p),0.,p) ; z=where(ismissing(z),0.,z) ; ps=where(ismissing(ps),0.,ps) ; sev=where(ismissing(sev),0.,sev) ; finally, accumulate values ; t_mam = t_mam + t ; q_mam = q_mam + q ; u_mam = u_mam + u ; v_mam = v_mam + v ; p_mam = p_mam + p ; ps_mam = ps_mam + ps ; z_mam = z_mam + z ; sev_mam = sev_mam + sev ; cape_mam = cape_mam + cape ; end if ; add the same variables here... if (nm.ge.5.and.nm.le.7) then ; same procedure for jja ; choose any of the variables to determine the nonmissing values ;XMISS::xmiss(t(0,:,:), miss_jja, nlon, nlat) XMISS::xmiss(t(0,:,:), miss_jja, nlon, nlat) ; now, go through and set missing values to zero t=where(ismissing(t),0.,t) q=where(ismissing(q),0.,q) u=where(ismissing(u),0.,u) v=where(ismissing(v),0.,v) ; p=where(ismissing(p),0.,p) z=where(ismissing(z),0.,z) ps=where(ismissing(ps),0.,ps) ; sev=where(ismissing(sev),0.,sev) t_jja = t_jja + t q_jja = q_jja + q u_jja = u_jja + u v_jja = v_jja + v ; p_jja = p_jja + p ps_jja = ps_jja + ps z_jja = z_jja + z ; sev_jja = sev_jja + sev cape_jja = cape_jja + cape end if ;copy_VarCoords(tt(0,0,:,:),sev_mam) ;copy_VarCoords(tt(:,:,:),t_mam) ;copy_VarCoords(tt(:,:,:),q_mam) ;copy_VarCoords(tt(:,:,:),u_mam) ;copy_VarCoords(tt(:,:,:),v_mam) ;copy_VarCoords(tt(:,:,:),p_mam) ;copy_VarCoords(tt(:,:,:),z_mam) ;copy_VarCoords(tt(0,:,:),miss_mam) ;copy_VarCoords(tt(0,:,:),ps_mam) ;copy_VarCoords(tt(0,0,:,:),sev_jja) copy_VarCoords(tt(:,:,:),t_jja) copy_VarCoords(tt(:,:,:),q_jja) copy_VarCoords(tt(:,:,:),u_jja) copy_VarCoords(tt(:,:,:),v_jja) ;copy_VarCoords(tt(:,:,:),p_jja) copy_VarCoords(tt(:,:,:),z_jja) copy_VarCoords(tt(0,:,:),miss_jja) copy_VarCoords(tt(0,:,:),ps_jja) ;sev_mam@units = "days" ;sev_mam@long_name = "accum occurrences" ; plot = gsn_csm_contour_map(wks,t_mam(nz,:,:), res) ; create plot ; plot = gsn_csm_contour_map(wks,sev_mam, res) ; create plot ; plot = gsn_csm_contour_map(wks,miss_mam, res) ; create plot else print("missing file="+file3) end if delete(f3) end do end do ; final averages ; to prevent division by zero, set the miss-counts to missing if equal to zero ;miss_mam = mask(miss_mam,miss_mam.eq.0., False) miss_jja = mask(miss_jja,miss_jja.eq.0., False) do nn=0,nlev-1 ;t_mam(nn,:,:) =t_mam(nn,:,:)/miss_mam ;q_mam(nn,:,:) =q_mam(nn,:,:)/miss_mam ;u_mam(nn,:,:) =u_mam(nn,:,:)/miss_mam ;v_mam(nn,:,:) =v_mam(nn,:,:)/miss_mam ;p_mam(nn,:,:) =p_mam(nn,:,:)/miss_mam ;z_mam(nn,:,:) =z_mam(nn,:,:)/miss_mam t_jja(nn,:,:) =t_jja(nn,:,:)/miss_jja q_jja(nn,:,:) =q_jja(nn,:,:)/miss_jja u_jja(nn,:,:) =u_jja(nn,:,:)/miss_jja v_jja(nn,:,:) =v_jja(nn,:,:)/miss_jja ;p_jja(nn,:,:) =p_jja(nn,:,:)/miss_jja z_jja(nn,:,:) =z_jja(nn,:,:)/miss_jja end do ;ps_mam =ps_mam/miss_mam ;cape_mam =cape_mam/miss_mam ps_jja =ps_jja/miss_jja cape_jja =cape_jja/miss_jja ; note: for the "sev", just accumulate, so that we have a total number of ; occurrences for the period (rather than an average) ; copy some metadata... copy_VarCoords(tt(0,:,:),cape_jja) ;copy_VarCoords(tt(:,:,0),cape_jja) cape_jja@units = "J/kg" cape_jja@long_name = "cape" ;copy_VarCoords(tt(0,0,:,:),cape_mam) ;cape_mam@units = "J/kg" ;cape_mam@long_name = "cape" ;copy_VarCoords(tt(0,0,:,:),sev_mam) ;sev_mam@units = "days" ;sev_mam@long_name = "occurrences" ;copy_VarCoords(tt(0,0,:,:),sev_jja) ;sev_jja@units = "days" ;sev_jja@long_name = "occurrences" copy_VarCoords(tt(:,:,:),z_jja) ;copy_VarCoords(tt(0,:,:,:),z_mam) ; dump out data ;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ; now, dump out the monthly averages to a ; netcdf file ;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ncfile = "jjafields.nc" system ("rm " + ncfile) ; remove any pre-existing file ncdf = addfile(ncfile ,"c") ; open output netCDF file ;=================================================================== ; create global attributes of the file ;=================================================================== fAtt = True ; assign file attributes fAtt@title = "Summer Averages" fAtt@domxmin = lon(0) fAtt@domxmax = lon(nlon-1) fAtt@domymin = lat(0) fAtt@domymax = lat(nlat-1) fileattdef( ncdf, fAtt ) ; define file attributes dimNames = (/"lat","lon"/) ; dimNames = (/"lon","lat"/) dimSizes = (/nlat,nlon/) ; dimSizes = (/nlon,nlat/) dimUnlim = (/False,False/) filedimdef(ncdf,dimNames,dimSizes,dimUnlim) ncdf->u_jja = u_jja ; lat, lon ncdf->v_jja = v_jja ; lat, lon ; ncdf->p_jja = p_jja ; lat, lon ncdf->ps_jja = ps_jja ; lat, lon ncdf->z_jja = z_jja ; lat, lon ncdf->t_jja = t_jja ; lat, lon ncdf->q_jja = q_jja ; lat, lon ncdf->lon = lon ncdf->lat = lat ; ncfile2 = "mamfields.nc" ; system ("rm " + ncfile2) ; remove any pre-existing file ; ncdf2 = addfile(ncfile2 ,"c") ; open output netCDF file ;=================================================================== ; create global attributes of the file ;=================================================================== ; fAtt = True ; assign file attributes ; fAtt@title = "Spring Averages" ; fAtt@domxmin = lon(0) ; fAtt@domxmax = lon(nlon-1) ; fAtt@domymin = lat(0) ; fAtt@domymax = lat(nlat-1) ; fileattdef( ncdf2, fAtt ) ; define file attributes ; dimNames = (/"lat","lon"/) ; dimSizes = (/nlat,nlon/) ; dimUnlim = (/False,False/) ; filedimdef(ncdf2,dimNames,dimSizes,dimUnlim) ; ncdf2->u_mam = u_mam ; lat, lon ; ncdf2->v_mam = v_mam ; lat, lon ; ncdf2->p_mam = p_mam ; lat, lon ; ncdf2->ps_mam = ps_mam ; lat, lon ; ncdf2->z_mam = z_mam ; lat, lon ; ncdf2->t_mam = t_mam ; lat, lon ; ncdf2->q_mam = q_mam ; lat, lon ;ncdf2->lon = lon ;ncdf2->lat = lat ; res@cnLevelSelectionMode = "ManualLevels" ; res@cnLevelSpacingF = 5. ; res@cnMinLevelValF = 0. ; res@cnMaxLevelValF = 50. ; plot = gsn_csm_contour_map(wks,cape_mam, res) ; create plot ; plot = gsn_csm_contour_map(wks,sev_mam, res) ; create plot ; plot = gsn_csm_contour_map(wks,u_mam(nz,:,:), res) ; create plot ; plot = gsn_csm_contour_map(wks,v_mam(nz,:,:), res) ; create plot ; plot = gsn_csm_contour_map(wks,t_mam(nz,:,:), res) ; create plot ; plot = gsn_csm_contour_map(wks,q_mam(nz,:,:), res) ; create plot ; plot = gsn_csm_contour_map(wks,p_mam(nz,:,:), res) ; create plot ; plot = gsn_csm_contour_map(wks,z_mam(nz,:,:), res) ; create plot ; plot = gsn_csm_contour_map(wks,ps_mam, res) ; create plot ; plot = gsn_csm_contour_map(wks,miss_mam, res) ; create plot plot = gsn_csm_contour_map(wks,cape_jja, res) ; create plot ; plot = gsn_csm_contour_map(wks,sev_jja, res) ; create plot plot = gsn_csm_contour_map(wks,u_jja(nz,:,:), res) ; create plot plot = gsn_csm_contour_map(wks,v_jja(nz,:,:), res) ; create plot plot = gsn_csm_contour_map(wks,t_jja(nz,:,:), res) ; create plot plot = gsn_csm_contour_map(wks,q_jja(nz,:,:), res) ; create plot ;plot = gsn_csm_contour_map(wks,p_jja(nz,:,:), res) ; create plot plot = gsn_csm_contour_map(wks,z_jja(nz,:,:), res) ; create plot plot = gsn_csm_contour_map(wks,ps_jja, res) ; create plot plot = gsn_csm_contour_map(wks,miss_jja, res) ; create plot end