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" ;***************************************************** ; ------------- MAIN DRIVER ------------------------ ;***************************************************** ;***************************************************** ; Miscellaneous initializations ;***************************************************** setfileoption("nc", "SuppressClose", False) yrStrt = 1990 yrLast = 2010 yrLast = yrStrt nmoStrt = 1 nmoLast = 12 var = "R_GDS4_ISBL" dRoot = "/data1/data/idosam/Era_netcdf2/" fRoot = "ei.oper.an.pl.regn128sc." ; file root print("===================================================================") print("dRoot="+dRoot) print("fRoot="+fRoot) print("===================================================================") ;***************************************************** ; Missing information from the file ;***************************************************** lev = (/1,2,3,5,7,10,20,30,50,70,100 \ ,125,150,175,200,225,250,300,350,400,450,500 \ ,550,600,650,700,750,775,800 ,825,850,875,900 \ ,925,950,975,1000 /) lev!0 = "lev" lev@long_name = "pressure level" lev@units = "hPa" nlat = 256 lat = latGau(nlat, "lat", "latitude", "degrees_north") mlon = 2*nlat lon = lonGlobeF(mlon, "lon", "longitude", "degrees_east") ;***************************************************** ; Read 4x-daily data ;***************************************************** do year=yrStrt,yrLast do nmo=nmoStrt,nmoLast yyyymm = year*100+nmo print(" ") print("=====> Calculating Monthly Mean: "+yyyymm+" <=====") diri = dRoot fili = systemfunc("cd "+diri+" ; ls "+fRoot+yyyymm+"*") print(fili) nfili = dimsizes(fili) f = addfiles( diri+fili, "r") ListSetType(f, "join") x = f[:]->$var$ printVarSummary(x) delete(x@initial_time) ; not nneded for monthly mean ; create time coordinate information date = toint(str_get_cols(fili, 24, 33)) yr = toint(str_get_cols(fili, 24, 27)) mon = toint(str_get_cols(fili, 28, 29)) day = toint(str_get_cols(fili, 30, 31)) hour = toint(str_get_cols(fili, 32, 33)) mn = new (dimsizes(hour), "integer") sc = new (dimsizes(hour), "integer") mn = 0 sc = 0 units = "hours since 1901-1-1 00:00:0.0" time = cd_inv_calendar(yr,mon,day,hour,mn,sc,units, 0) time!0= "time" time&time = time delete(time@_FillValue) x!0 = "time" ; create time coordinate x&time = time x!1 = "lev" ; rename to nicer dimension names x!2 = "lat" x!3 = "lon" xMon = dim_avg_n_Wrap(x, 0) printVarSummary(xMon) ; [lev | 37] x [lat | 256] x [lon | 512] ; add a time dimension dimx = dimsizes(x) rankx = dimsizes(dimx) print(dimx) print("rank="+rankx) dimMon = dimx dimMon(0) = 1 print(dimMon) if (rankx.eq.3) then xMON = conform_dims( dimMon, xMon, (/1,2/) ) copy_VarMeta(xMon, xMON(0,:,:)) else xMON = conform_dims( dimMon, xMon, (/1,2,3/) ) copy_VarMeta(xMon, xMON(0,:,:,:)) end if ;;printVarSummary(xMON) xMON!0 = "time" xMON&time = x&$x!0$(0) ; assign is day of month xMON!1 = "lev" xMON!2 = "lat" xMON!3 = "lon" xMON&lev = lev xMON&lat = lat xMON&lon = lon printVarSummary(xMON) yyyymm = date(0)/10000 yyyymm!0 = "time" yyyymm@units= "yyyymm" printVarSummary(yyyymm) ;***************************************************** ; Write monthly means to netCDF ; pleas read: http://www.ncl.ucar.edu/Applications/method_1.shtml ;***************************************************** fout = "RHMON."+date(0)+".nc" print("Writing out netcdf: "+fout) system("/bin/rm -f "+fout) ncdf = addfile(fout, "c") filedimdef(ncdf,"time",-1, True) ncdf->date = yyyymm ; only the ncdf->RH = xMON delete([/ fili, x, xMon, xMON, dimx, dimMon /]) delete([/ yr, mon, day, hour, mn, sc, date, time /]) end do ; nmo end do ; year