Hi,
I am trying to calculate monthly averages for some variables and I am
still getting this error and don't know why.
fatal:NclMalloc Failed:[errno=12]
fatal:New: could not create new array:[errno=12]
fatal:Execute: Error occurred at or near line 193 in file simple_new.ncl
My script is the following:
; in this script, we compute CAPE, shear, etc., and then evaluate
; the frequency at which threshholds of these parameters are met
; this version is relevant for the CAM runs
; this particular version takes advantage of the NCL command line arguments
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 SHEAR "./newshear.so"
;external DISCRIM "./newdiscrim.so"
;external DISCRIM "./newdiscrim_surface.so"
external DISCRIM "./newdiscrim2.so"
; have this now set up to run through a decade at a time
; just need to change the first line, and also the paths, as needed
begin
thresh_shr = 0.
thresh_cape = 1000.
; comment out these two lines when running this in batch mode
; (the decade and filstrng will be provided at the command line in the
; shell script)
decade = 1991
filstrng = "030bES_T85"
; path to the binary files
; trying out new scratch space
main1 = "/scratch/scratch96/c/cvillanu/"
main3 = "/scratch/scratch96/c/cvillanu/"
; for the output files...
monstrng =
(/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"/)
file3ext = (/"-01-01-","-03-15-","-05-27-","-08-08-","-10-20-"/)
; the number of days per month ... non leap year
jday = (/0,31,59,90,120,151,181,212,243,273,304,334,365/)
; constants
grav = 9.80616
pt = 50 ; in hPa
;************************************************
; create plot
;************************************************
; wks = gsn_open_wks("pdf","color") ; open a ps file
wks = gsn_open_wks("ncgm","color") ; open a ps file
; wks = gsn_open_wks("x11","color") ; open a ps file
gsn_define_colormap (wks,"gui_default") ; choose color map
;******************************************************
; get some initial information from an arbitrary netCDF file
;******************************************************
;Pointers
; the cam2.h0 are the monthly averaged files, but contain lots other fields
file1=filstrng+".cam2.h1."+decade+"-01-01-00000.nc"
file0=filstrng+".cam2.h0."+decade+"-01.nc"
file2=filstrng+".clm2.h0."+decade+"-01.nc"
f = addfile(main1+file1, "r")
f2 = addfile(main1+file2, "r")
f0 = addfile(main1+file0, "r")
print("file1="+file1)
;Read Variables
l1 = f->lon
lon=doubletofloat(l1) ; coerce, for future calculations
l2 = f->lat
lat=doubletofloat(l2) ; coerce, for future calculations
p0 = f->P0 ; reference pressure (Pa)
P0=doubletofloat(p0) ; coerce, for future calculations
;sanity check...
tt = f->T(:,4:25,:,:)
print(dimsizes(tt))
;print(dims)
lev = f->lev(4:25) ; lev is hybrid sigma levelts (at midpoints)
;lev = f->lev ; lev is hybrid sigma levelts (at midpoints)
HT = f2->ZBOT(0,:,:) ; Terrain Height (lat, lon)
LU = f2->landmask ; land/ocean mask (0.=ocean and 1.=land)
; try this instead...
PHIS = f0->PHIS(0,:,:) ; surface geopotential (m2/s2)
; (monthly average)
HT = PHIS/grav ; proxy for height of topography
; a few others...
a1 = f->hyam(4:25) ; hybrid A coefficient at layer
midpoints
;a1 = f->hyam ; hybrid A coefficient at layer midpoints
hyam=doubletofloat(a1) ; coerce, for future calculations
b1 = f->hybm(4:25) ; hybrid B coefficient at layer midpoints
;b1 = f->hybm ; hybrid B coefficient at layer midpoints
hybm=doubletofloat(b1) ; coerce, for future calculations
; this is done to eliminate the missing value flag of zero, which
; otherwise sets subquent calculations to zero
delete(HT@_FillValue)
nlat = dimsizes(lat)
nlon = dimsizes(lon)
print("nlat="+nlat)
print("nlon="+nlon)
nlev = dimsizes(lev)
nz=nlev-1
print("nlev="+nlev)
; assign these attributes
cape = new((/nlat,nlon/), float)
ripcape = new((/4,nlat,nlon/), float)
s06 = new((/nlat,nlon/), float)
uavg = new((/nlev,nlat,nlon/), float)
vavg = new((/nlev,nlat,nlon/), float)
tavg = new((/nlev,nlat,nlon/), float)
qavg = new((/nlev,nlat,nlon/), float)
zavg = new((/nlev,nlat,nlon/), float)
pavg = new((/nlev,nlat,nlon/), float)
ppsavg = new((/nlat,nlon/), float) ; surface pressure, get surface z and
q
;tsavg = new((/nlat,nlon/), float)
;qsavg = new((/nlat,nlon/), float)
sev = new((/nlat,nlon/), float)
copy_VarCoords(tt(0,:,:,:),uavg)
uavg_at_units = "m/s"
uavg_at_long_name = "mean zonal velocity"
copy_VarCoords(tt(0,:,:,:),vavg)
vavg_at_units = "m/s"
vavg_at_long_name = "mean meridional velocity"
copy_VarCoords(tt(0,:,:,:),tavg)
tavg_at_units = "K"
tavg_at_long_name = "mean temperature"
copy_VarCoords(tt(0,:,:,:),qavg)
qavg_at_units = "kg/kg"
qavg_at_long_name = "mean specific humidity"
copy_VarCoords(tt(0,:,:,:),zavg)
zavg_at_units = "m"
zavg_at_long_name = "mean geopotential height"
copy_VarCoords(tt(0,:,:,:),pavg)
pavg_at_units = "Pa"
pavg_at_long_name = "mean pressure"
copy_VarCoords(tt(0,0,:,:),ppsavg)
ppsavg_at_units = "Pa"
ppsavg_at_long_name = "Surface Pressure"
;copy_VarCoords(tt(0,0,:,:),tsavg)
;tsavg_at_units = "K"
;tsavg_at_long_name = "2m mean air temperature"
;copy_VarCoords(tt(0,0,:,:),qsavg)
;qsavg_at_units = "kg/kg"
;qsavg_at_long_name = "2m mean specific humidity"
copy_VarCoords(tt(0,0,:,:),sev)
sev_at_units = "days"
sev_at_long_name = "SEV"
copy_VarCoords(tt(0,0,:,:),cape)
cape_at_units = "J/kg"
cape_at_long_name = "CAPE"
copy_VarCoords(tt(0,0,:,:),ripcape(0,:,:))
ripcape_at_units = "J/kg"
ripcape_at_long_name = "CAPE"
; try defining the variables
T = new((/365,nlev,nlat,nlon/), float)
Q = new((/365,nlev,nlat,nlon/), float)
U = new((/365,nlev,nlat,nlon/), float)
V = new((/365,nlev,nlat,nlon/), float)
Z3 = new((/365,nlev,nlat,nlon/), float)
pps = new((/365,nlat,nlon/), float)
rc = new((/365,nlat,nlon/), float)
;tts = new((/365,nlat,nlon/), float)
;qqs = new((/365,nlat,nlon/), float)
; create 3D array of pressure
P = new((/nlev,nlat,nlon/), float)
ntot = 1460 ; this is the sum of all 6-hr periods in a 365-day year
; fortunately, this run assumes no leap years
; in this version, only do a year at a time
do nd = 0,0
year = decade+nd
tbeg = 0
tend = 72
;do nf = 0,4
; only read the files through Aug
do nf = 0,3
file3=main3+filstrng+".cam2.h1."+year+file3ext(nf)+"00000.nc"
file4=main3+filstrng+".cam2.h2."+year+file3ext(nf)+"00000.nc"
fa = addfile(file3, "r")
fb = addfile(file4, "r")
print("file3="+file3)
print("tbeg,tend="+tbeg+" "+tend)
;Read Variables
T(tbeg:tend,:,:,:) = fa->T(::4,4:25,:,:) ; here, the ::4 indicates
; every 4th time, i.e., 00 UTC
Q(tbeg:tend,:,:,:) = fa->Q(::4,4:25,:,:)
U(tbeg:tend,:,:,:) = fa->U(::4,4:25,:,:)
V(tbeg:tend,:,:,:) = fa->V(::4,4:25,:,:)
Z3(tbeg:tend,:,:,:) = fa->Z3(::4,4:25,:,:)
pps(tbeg:tend,:,:) = fa->PS(::4,:,:)
; Note** need to read the new variables!
; for example,
;qqs(tbeg:tend,:,:) = f2->Q2M(0,:,:)
;tts(tbeg:tend,:,:) = f2->TSA(0,:,:)
delete(fa)
delete(fb)
tbeg=tend+1
tend = tbeg+72
end do
; start processing the months...
; only do the warm season months MAM - JJA
do nm = 6, 6
;do nm = 2, 7
print("the month is"+nm)
; initialize to zero for each *month*
uavg = 0.
vavg = 0.
tavg = 0.
qavg = 0.
zavg = 0.
pavg = 0.
ppsavg = 0.
;tsavg = 0
;qsavg = 0
sev = 0
ntlev = 0
dybeg = jday(nm)
dyend = jday(nm+1)-1
print("for month="+nm+" dybeg="+dybeg+" dyend="+dyend)
do nt = dybeg,dyend ; this should read only the 00 UTC files
T3D = T(nt,:,:,:) ; as currently configured, that's the only
Q3D = Q(nt,:,:,:) ; time in the arrays
U3D = U(nt,:,:,:)
V3D = V(nt,:,:,:)
ps = pps(nt,:,:)
GP = Z3(nt,:,:,:)
;ts = tts(nt,:,:)
;qs = qqs(nt,:,:)
do k=0,nz
P(k,:,:)=0.01*(hyam(k)*P0+hybm(k)*ps)
end do
ps=ps*.01
ntlev=ntlev+1
; now to the CAPE
; for cam, arrays are already top to bottom
;ripcape = wrf_cape_2d(P,T3D,Q3D,GP,HT,ps,True)
ripcape = rip_cape_2d(P,T3D,Q3D,GP,HT,ps,True) ; the True flag indicates
; sigma level data
; since this is a hybrid
coord
; will need to be careful...
cape = ripcape(0,:,:)
; let's try plotting some of the variables...
;************************************************
; create plot
;************************************************
res = True ; plot mods desired
res_at_gsnSpreadColors = True
res_at_cnFillOn = True ; turn on color fill
res_at_mpFillOn = False ; turn off continent gray
res_at_mpOutlineBoundarySets = "GeophysicalAndUSStates"
res_at_cnLinesOn = True
res_at_cnLineLabelsOn = True
res_at_gsnAddCyclic = False ; this applies to the windowed
domain
res_at_mpLimitMode = "LatLon" ; control map zoom with lat/lons
res_at_mpMaxLatF = 55. ; actual values of zoom
res_at_mpMinLatF = 25.
res_at_mpMinLonF = 220.
res_at_mpMaxLonF = 300.
res_at_gsnMaximize = False
plot = gsn_csm_contour_map(wks,cape, res) ; create plot
;****************************************************************
; New: compute 0-6km shear, using external fortran function
;****************************************************************
SHEAR::bshear(U3D, V3D, GP, HT, s06, nlon, nlat, nlev)
; find the environment days (i.e., those that exceed the shear & cape
thresh)
; include the ppsavg, and avg sfc variables in this call...
;Note**: you need to pass the original arrays of ps, ts, qs to discrim
; and then you need to modify discrim to accomodate this
DISCRIM::discrim(cape,s06,U3D,V3D,T3D,Q3D,GP,P,sev,uavg,vavg,tavg,qavg,zavg,pavg,ppsavg,nlon,nlat,nlev,thresh_shr,thresh_cape)
end do ; all days in a month
print("getting averages...ntlev="+ntlev)
; before computing final avg...set avg to missing where sev = 0
sev = mask(sev,sev.eq.0, False)
do k = 0,nz
uavg(k,:,:) = uavg(k,:,:)/sev(:,:)
vavg(k,:,:) = vavg(k,:,:)/sev(:,:)
tavg(k,:,:) = tavg(k,:,:)/sev(:,:)
qavg(k,:,:) = qavg(k,:,:)/sev(:,:)
pavg(k,:,:) = pavg(k,:,:)/sev(:,:)
zavg(k,:,:) = zavg(k,:,:)/sev(:,:)
ppsavg(:,:) = ppsavg(:,:)/sev(:,:)
; tsavg(:,:) = tsavg(:,:)/sev(:,:)
; qsavg(:,:) = qsavg(:,:)/sev(:,:)
; include ppsavg and others here...
end do
; note: when you begin your 'production' processing, turn off plotting
plot = gsn_csm_contour_map(wks,sev, res) ; create plot
plot = gsn_csm_contour_map(wks,ppsavg, res) ; create plot
; plot = gsn_csm_contour_map(wks,tsavg, res) ; create plot
; plot = gsn_csm_contour_map(wks,qsavg, res) ; create plot
; plot = gsn_csm_contour_map(wks,vavg(nz,:,:), res) ; create plot
; plot = gsn_csm_contour_map(wks,zavg(nz,:,:), res) ; create plot
; plot = gsn_csm_contour_map(wks,pavg(nz,:,:), res) ; create plot
; plot average cape
; plot = gsn_csm_contour_map(wks,avgcape, res) ; create plot
; plot = gsn_csm_contour_map(wks,mavgcape, res) ; create plot
; plot = gsn_csm_contour_map(wks,sev, res) ; create plot
;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
; now, dump out the monthly averages to a
; netcdf file
;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ncfile = year + monstrng(nm) + ".nc"
ncdf = addfile(ncfile ,"c") ; open output netCDF file
;===================================================================
; create global attributes of the file
;===================================================================
fAtt = True ; assign file attributes
fAtt_at_title = "CAM output"
fAtt_at_domxmin = lon(0)
fAtt_at_domxmax = lon(nlon-1)
fAtt_at_domymin = lat(0)
fAtt_at_domymax = lat(nlat-1)
fileattdef( ncdf, fAtt ) ; define file attributes
dimNames = (/"lat","lon"/)
dimSizes = (/nlat,nlon/)
dimUnlim = (/False,False/)
filedimdef(ncdf,dimNames,dimSizes,dimUnlim)
ncdf->uavg = uavg ; nlev, lat, lon
ncdf->vavg = vavg ; nlev, lat, lon
ncdf->tavg = tavg ; nlev, lat, lon
ncdf->qavg = qavg ; nlev, lat, lon
ncdf->pavg = pavg ; nlev, lat, lon
ncdf->zavg = zavg ; nlev, lat, lon
ncdf->sev = sev
ncdf->ppsavg = ppsavg
;ncdf->tsavg = tsavg
;ncdf->qsavg = qsavg
; add other variables to the netcdf output
ncdf->lon = lon
ncdf->lat = lat
end do ; all months
end do ; all years
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Aug 02 2009 - 10:12:26 MDT
This archive was generated by hypermail 2.2.0 : Wed Aug 05 2009 - 20:36:03 MDT