Hello folks,
this will be my first contribution to ncl-talk, so please don't hesitate
to tell me what I did not include in my message ;-).
I am processing satellite level2 data in order to put it on a normal
global grid and I am using code that was provided to me by Dennis. The
satellite data normally has array sizes of 4224*191*float and the output
array size is 721*1441*float. After many iterations, the error message:
fatal:NclMalloc Failed:[errno=12]
Segmentation fault
appears (out of the blue for me ;-) ). From reading other ncl-talk
messages (this error seems to happen quite often), it is pretty sure
that this is a memory issue. I am pretty much deleting all variables
created for every iteration in order for them to not grow over the 2GB
limit. I will attach my code to this message (the principal code starts
at approx. line 107).
I hope you guys can help me handle these memory issues ! And by the way,
is there some other way to handle these NCL memory problems besides
deleting variables all the time and configurating the memory capacities
of the machine NCL is running on ?
Regards and many thanks,
Karsten
-- Karsten Peters Meteorological Institute University of Hamburg c/o Max-Planck Institute for Meteorology Bundesstrasse 53 D-20146 Hamburg ------------------------------------------------------------------ Tel.: ++49 (0)40 41173 147 Fax.: ++49 (0)40 41173 298 Email: karsten.peters_at_zmaw.de URL: http://www.mpimet.mpg.de/~peters.karsten
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"
external BINF "./fbindata_atsr.so"
;---
function parseFileName ( fNam:string ) ;filename(s) given to function
; parse the file names and extract data information
; 1 2
; 0123456789012345678901234567890
; GRAPE_sepVarsl2_200001092200.nc
local onNam, n, output, tmp_c
begin
nNam = dimsizes( fNam )
if (nNam.eq.1) then
output = new( 5,integer,"No_FillValue")
tmp_c = stringtochar(fNam)
output(0) = stringtointeger((/tmp_c(16:19)/)) ; YYYY
output(1) = stringtointeger((/tmp_c(20:21)/)) ; MM
output(2) = stringtointeger((/tmp_c(22:23)/)) ; DD
output(3) = stringtointeger((/tmp_c(24:25)/)) ; HH
output(4) = stringtointeger((/tmp_c(26:27)/)) ; MN
else ;array with date information for every file
output = new( 5,integer,"No_FillValue")
do n=0,nNam-1
tmp_c = stringtochar(fNam(n))
output(0) = stringtointeger((/tmp_c(16:19)/)) ; YYYY
output(1) = stringtointeger((/tmp_c(20:21)/)) ; MM
output(2) = stringtointeger((/tmp_c(22:23)/)) ; DD
output(3) = stringtointeger((/tmp_c(24:25)/)) ; HH
output(4) = stringtointeger((/tmp_c(26:27)/)) ; MN
end do
end if
return (output)
end
;---
function create3d (x[*][*]:numeric, time[*]:numeric)
; convert a 2-dimensional to a 3-dimensional variable
; with time as a coordinate variable.
local dimx, ntim, nlat, mlon, x3d
begin
dimx = dimsizes(x)
ntim = dimsizes(time)
nlat = dimx(0)
mlon = dimx(1)
; perform expansion
x3d = conform_dims( (/ntim,nlat,mlon/), x, (/1,2/)) ;nlat and nlon correspond to dimension index 1 and 2 on x
copy_VarAtts(x, x3d)
x3d!0 = "time"
x3d&time= time
copy_VarCoords(x, x3d(0,:,:) ) ;copies all named dimensions and coordinate variables
return(x3d)
end
;++++++++++++++++++++++++ MAIN +++++++++++++++++++++++++++
begin
;************************************************************
;create the variable array
;************************************************************
vars = (/"Sol_zen_angle", "Sat_ zen_angle", "Rel_az_angle","Refl_var_055", "Refl_var_067", "Refl_var_087", "Refl_var_16", "Refl_var_37", "BT_var_11", "BT_var_12","Refl_055", "Refl_067", "Refl_087", "Refl_16", "Refl_37", "BT_11", "BT_12","log_AOD", "log_aer_eff_rad", "Surface_Alb_055","log_AOD_uncty", "log_aer_eff_rad_uncty", "Surface_Alb_055_uncty","Surface_Alb_067", "Surface_Alb_087", "Surface_Alb_16","log_COD", "log_c_eff_rad", "CTP", "CF", "ST", "CTH", "CTT", "LWP","log_COD_uncty", "log_c_eff_rad_uncty", "CTP_uncty", "CF_uncty", "ST_uncty", "CTH_uncty", "CTT_uncty", "LWP_uncty"/)
vNums = dimsizes(vars)
PLOT = True
netCDF= True
;*****************************************************************
; grid definition
;*****************************************************************
nlatb = 721 ; grid *boundaries*
mlonb = 1441
latb = fspan( -90, 90, nlatb) ; -90, 60
lonb = fspan(-180, 180, mlonb)
dlon = lonb(1)-lonb(0) ; grid spacing
dlat = latb(1)-latb(0)
nlat = nlatb-1 ; grid coordinates
mlon = mlonb-1
lat = latGlobeFo(nlat,"lat","latitude","degrees_north")
lon = lonGlobeFo(mlon,"lon","longitude","degrees_east")
lon = (/ lon - 180. /) ; subtract 180 from all values
lon&lon = lon ; update coordinates
;*****************************************************************
; Variables to hold binned quantities
;*****************************************************************
GBIN = new ( (/nlat,mlon/), float, -999 ) ;corresponds to my histogram
GKNT = new ( (/nlat,mlon/), float, -999 ) ;correspons to my histogram_count
GBIN = 0.0 ; initialization
GKNT = 0.0
PLOT = False
netCDF= True
byear = 2000
eyear = 2001
;do yr=byear,eyear
yr=2000
cyr=sprinti("%0.4i",yr)
bmon = 1
emon = 12
;do mn=bmon,emon
mn=1
cmn=sprinti("%0.2i",mn)
bday = 1
eday = 31
;do dd=bday,eday
dd=09
cday=sprinti("%0.2i",dd)
;*****************************************************************
; File info
;*****************************************************************
diri = "/scratch/local1/u231008/PhD/atsr_v2/"
c=" GRAPE_sepVarsl2_"+cyr+cmn+cday+"*.nc"
fili = systemfunc("cd "+diri+" ; ls"+c)
nfil = dimsizes( fili )
print("nfil="+nfil)
;*****************************************************************
; time/date
;*****************************************************************
ydhm = parseFileName( fili(0) )
yyyy = ydhm(0)
mn = ydhm(1)
day = ydhm(2)
hh = ydhm(3)
mi = ydhm(4)
;********************************************************************
;*****************************************************************
; netCDF
;*****************************************************************
if (netCDF) then ;in this case, netCDF = True
sec = 0.0
; COARDS/CF time
tunits = "hours since 1990-1-1 00:00:0.0"
time = ut_inv_calendar(yyyy,mn,day,hh,mi,sec,tunits, 0)
time!0 = "time"
time_at_units = tunits
; gregorian date
date = yyyy*10000 + mn*100 + day
date!0 = "time"
date_at_units = "yyyymmdd"
diro = "./" ; output directory
filo = "GRAPE_atsr-2_quart_"+date ; output file name
fout = diro+filo+".nc"
system("/bin/rm -f "+fout) ; remove any pre-existing file
ncdf = addfile(fout ,"c") ; open output netCDF file
;================================================================
; create global attributes of the file (not required)
;================================================================
fAtt = True ; assign file attributes
fAtt_at_title = "GRAPE ATSR-2 Data gridded to 0.25x0.25: " + yyyy+mn+cday
fAtt_at_source_file = "GRAPE_sepVarsl2_"+yyyy+mn+day+"hhmn.nc"
;fAtt_at_Conventions = "CF 1.0"
fAtt_at_creation_date = systemfunc ("date")
fileattdef( ncdf, fAtt ) ; copy file attributes
; recommended (generally)
filedimdef(ncdf,"time",-1,True) ; make time UNLIMITED dimension
ncdf->time = time
ncdf->date = date
end if
;**********************************************************************
tStrt = systemfunc("date") ; time the loop (wall clock)
do k=0,vNums-1
;*****************************************************************
; Variables to hold binned quantities
;*****************************************************************
GBIN = new ( (/nlat,mlon/), float) ;corresponds to my histogram
GKNT = new ( (/nlat,mlon/), float) ;correspons to my histogram_count
GBIN = 0.0 ; initialization
GKNT = 0.0
print(""+vars(k))
do nf=0,nfil-1
print(nf+" "+fili(nf))
f = addfile(diri+fili(nf), "r")
; read data
x = f->$vars(k)$
xdims = dimsizes(x)
if (nf .eq. 0) then
xMeta = new((/xdims(0),xdims(1)/),float,-999)
copy_VarMeta(x,xMeta)
end if
lat2d = f->Lat
lon2d = f->Lon
; reshape to use "ind"
x1d = ndtooned( x )
lat1d = ndtooned( lat2d )
lon1d = ndtooned( lon2d )
; good (non-missing) values
iGood = ind(.not.ismissing(x1d))
if (.not.ismissing(iGood(0))) then
z = x1d(iGood)
zlat = lat1d(iGood)
zlon = lon1d(iGood)
nz = dimsizes( z )
; grid subscript indices
n = numeric2int( (zlat-latb(0))/dlat, 0) ;one dim. field of truncated integer indices
m = numeric2int( (zlon-lonb(0))/dlon, 0) ;of dimension nz
;Fortran for speed
BINF::bindata_atsr(nlat,mlon,GBIN,GKNT, nz,z, n,m)
delete( z ) ; size may change
delete( zlat)
delete( zlon)
delete( n )
delete( m )
delete ( nz)
end if
delete( x ) ; size may change
delete(lat2d)
delete(lon2d)
delete( x1d)
delete(lat1d)
delete(lon1d)
delete(iGood)
delete(f)
end do ;files
;*****************************************************************
; Perform averaging
;*****************************************************************
; avoid division by 0.0
GKNT = where(GKNT.eq.0 , GKNT@_FillValue, GKNT)
GBIN = GBIN/GKNT
copy_VarAtts(xMeta,GBIN)
delete(xMeta)
GBIN!0 = "lat"
GBIN!1 = "lon"
GBIN&lat = lat
GBIN&lon = lon
copy_VarCoords(GBIN, GKNT) ; copy coords
if (isatt(GBIN,"long_name")) then
ln = GBIN_at_long_name
GBIN_at_long_name = "BINNED: "+ln
GKNT_at_long_name = "BINNED COUNT: "+ln
end if
;if (isfilevaratt(f, vNam, "units")) then
; GBIN_at_units = f->$vNam$@units ;that is very very interesting !!!
;end if
;*****************************************************************
; time/date
;*****************************************************************
ydhm = parseFileName( fili(0) )
yyyy = ydhm(0)
ddd = ydhm(1)
hh = ydhm(2)
mn = ydhm(3)
yyyyddd = yyyy*1000 + ddd
monday = monthday(yyyy, ddd)
month = monday/100
day = monday - month*100 ; day of month
;*****************************************************************
; write netCDF
;*****************************************************************
if (netCDF) then
x3d = create3d( GBIN, time)
ncdf->$vars(k)$ = x3d
kNam = vars(k)+"_knt"
x3d = create3d( GKNT, time)
ncdf->$kNam$ = x3d
end if
delete(GBIN)
delete(GKNT)
end do ;variables
wallClockElapseTime(tStrt, "Main File Loop", 0)
;end do ;days
;end do ;months
;end do ;years
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Jan 14 2009 - 09:15:32 MST
This archive was generated by hypermail 2.2.0 : Fri Jan 16 2009 - 14:05:56 MST