[ncl-talk] problem in writing binary file ...
Ufuk Utku Turuncoglu (BE)
u.utku.turuncoglu at be.itu.edu.tr
Thu Sep 11 03:03:29 MDT 2014
Hi Dennis,
Thanks for your help. I could write the data correctly after making some
modifications in the fortran code (it is attached). In this case, i have
to set the binary file as "formatted" plus "sequential" to create file.
In this case, the fortran code read the file (generated by NCL) without
any problem. I don't know why but the fortran code confuses when it try
to read the binary file which is written in "unformatted" option. So, it
is still unclear for me. Anyway, my problem is solved. Thanks again for
your help.
Regards,
--ufuk
On 11/09/14 06:22, Dennis Shea wrote:
> NCL's binary write functions do not have the ability to write mixed
> type records.
> Further, fortran default sequential modepre/post-pend each record
> with a 'hidden' record separator .
>
> I'd suggest creating a fortran subroutine to do what you want.
> Untested f77
>
> C NCLFORTSTART
> subroutine wrhead(fname, llc, urc, nlat, mlon, dx, dy, kk)
> implicit none
> C input
> integer nlat, mlon, llc(2), kk
> real urc(2), dx, dy
> character*(*) fname
> C NCLEND
>
>
> open(unit=11, file=fname, form="unformatted", access="sequential")
>
> write(11) llc(0), urc(0), llc(1), urc(1), dx, dy, nlat, nlon, kk)
> return
> end
>
> C NCLFORTSTART
> subroutine wrrec(fname, ......)
> implicit none
> .....
> C NCLEND
>
> return
> end
>
>
> On Wed, Sep 10, 2014 at 5:54 PM, Jatin Kala <jatin.kala.jk at gmail.com
> <mailto:jatin.kala.jk at gmail.com>> wrote:
>
> Hi,
> This may or may not help solve your problem.
> My general experience is, it's easier to deal with direct-access
> binary
> versus sequential.
> I would consider using fbindirwrite rather than fbinrecwrite.
>
> Also in our OPEN statement, have you tried to specify the
> "ACCESS=" option?
> Cheers,
> Jatin.
>
>
> On 10/09/14 8:33 PM, Ufuk Utku Turuncoglu (BE) wrote:
> > Hi,
> >
> > I am trying to write a binary file, which is in following form,
> >
> > header 1 - 6 float and 2 integer
> > header 2 - character (len = 14), holds time information
> > data 1 - 2 dimensional data written as a chunk of lons (inside a
> lat loop)
> > data 2 - 2 dimensional data written as a chunk of lons (inside a
> lat loop)
> > header 2
> > data 1
> > data 2
> > ...
> > ...
> >
> > So, my ncl script writes the data as follows,
> >
> > ihead = (/ llc(0), urc(0), llc(1), urc(1), 0.25, 0.25, nlat,
> nlon, 3 /)
> > fbinrecwrite(fname, -1, ihead)
> >
> > do k = 0, ntime-1
> > print(k+" "+date_str(k))
> > fbinrecwrite(fname, -1, date_str(k))
> > do i = 0, nlat-1
> > fbinrecwrite(fname, -1, uwnd_rmap(k,i,:))
> > end do
> > do i = 0, nlat-1
> > fbinrecwrite(fname, -1, vwnd_rmap(k,i,:))
> > end do
> > end do
> >
> > The problem is that the fortran code try to read "header 1" with a
> > single read command. Like following,
> >
> > OPEN
> >
> (UNIT=IU01,FILE=FILE01(1:LEN),FORM='UNFORMATTED',STATUS='OLD',IOSTAT=IOS)
> > READ (IU01,IOSTAT=IOS) SOUTH, NORTH, WEST, EAST, D_LON,
> D_LAT,
> > N_LAT, N_LON
> >
> > but the data type are different for header (real for SOUTH,
> NORTH, WEST,
> > EAST, D_LON, D_LAT and integer for N_LAT, N_LON). As it can be seen
> > above, i tried to write the header as,
> >
> > ihead = (/ llc(0), urc(0), llc(1), urc(1), 0.25, 0.25, nlat,
> nlon, 3 /)
> >
> > but fortran file can not read it correctly. It reads "header 1"
> > correctly for SOUTH, NORTH, WEST, EAST, D_LON, D_LAT variables but
> > N_LAT, N_LON is wrong. So, i think that i have to find another
> way to
> > write header in a single record. I also check that NCL float and
> integer
> > data types and both of them are 32 bits but even the size are same
> > fortran reads them wrongly.
> >
> > I just wonder that is there any other way to write this kind of
> header
> > (has multiple data types) to the binary file.
> >
> > Regards,
> >
> > --ufuk
> >
> > _______________________________________________
> > ncl-talk mailing list
> > List instructions, subscriber options, unsubscribe:
> > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20140911/33f3f50a/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/applefile
Size: 433 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20140911/33f3f50a/attachment.bin
-------------- next part --------------
A non-text attachment was scrubbed...
Name: wam_utils.f
Type: application/octet-stream
Size: 2593 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20140911/33f3f50a/attachment.obj
-------------- next part --------------
;-----------------------------------------------------------
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
external WAM "./wam_utils.so"
;-----------------------------------------------------------
;***********************************************************
;*** Ufuk Turuncoglu ***************************************
;*** Send bug reports to turuncu at be.itu.edu.tr *************
;***********************************************************
begin
;---------------------------------------------------------
; User defined parameters
;---------------------------------------------------------
;--- input file for forcing ---
atm_in_file = "MED50_WIND_r46.nc"
;--- create grid definitions ---
fresh = True
;--- destination grid properties --
llc = (/ 29.0, -7.0 /)
urc = (/ 47.0, 38.0 /)
dx = "0.25deg"
;--- debug ---
debug = True
;---------------------------------------------------------
; ********************************************************
; *** Main section ***************************************
; ********************************************************
;
; Do not change the code after this point if it is not
; really neccecary !!!
;---------------------------------------------------------
;--- set file option ---
setfileoption("nc", "Format", "LargeFile")
;--- open atm model file ---
atm_nc = addfile(atm_in_file, "r")
atm_lat2d = atm_nc->xlat
atm_lon2d = atm_nc->xlon
atm_msk2d = atm_nc->mask
;--- get time dimension ---
time = atm_nc->time
date = tointeger(cd_calendar(time, 0))
date_str = sprinti("%04d", date(:,0))+\
sprinti("%02d", date(:,1))+\
sprinti("%02d", date(:,2))+\
sprinti("%02d", date(:,3))+\
sprinti("%02d", date(:,4))+\
sprinti("%02d", date(:,5))
ntime = dimsizes(time)
;---------------------------------------------------------
; STEP (1)
; Rotate wind
;---------------------------------------------------------
uwnd = atm_nc->uas(:,0,:,:)
vwnd = atm_nc->vas(:,0,:,:)
cone = tofloat(atm_nc at grid_factor)
clon = tofloat(atm_nc at longitude_of_projection_origin)
dumm = wrf_uvmet(uwnd, vwnd, atm_lat2d, atm_lon2d, clon, cone)
uwnd_rot = dumm(0,:,:,:)
vwnd_rot = dumm(1,:,:,:)
;---------------------------------------------------------
; STEP (2)
; Interpolate to regular lat-lon grid
;---------------------------------------------------------
;--- create SCRIP definition of source grid ---
sfile = "src_grd.nc"
if (.not. isfilepresent(sfile) .or. fresh) then
opt = True
opt at ForceOverwrite = True
opt at PrintTimings = True
opt at Title = "RCM 50km grid"
curvilinear_to_SCRIP(sfile, atm_lat2d, atm_lon2d, opt)
delete(opt)
end if
dfile = "dst_grd.nc"
if (.not. isfilepresent(dfile) .or. fresh) then
opt = True
opt at ForceOverwrite = True
opt at PrintTimings = True
opt at Title = "Regular grid"
opt at LLCorner = llc
opt at URCorner = urc
latlon_to_SCRIP(dfile, dx, opt)
delete(opt)
end if
;--- create weight matrix file ---
wfile = "wgt_src_to_dst.nc"
if (.not. isfilepresent(wfile) .or. fresh) then
opt = True
opt at ForceOverwrite = True
opt at PrintTimings = True
opt at SrcRegional = True
opt at DstRegional = True
opt at InterpMethod = "bilinear"
ESMF_regrid_gen_weights(sfile, dfile, wfile, opt)
end if
;--- perfrom interpolation ---
opt = True
opt at PrintTimings = True
uwnd_rmap = ESMF_regrid_with_weights(uwnd_rot, wfile, opt)
vwnd_rmap = ESMF_regrid_with_weights(vwnd_rot, wfile, opt)
lat1d = uwnd_rmap&lat
lon1d = uwnd_rmap&lon
dims = dimsizes(uwnd_rmap)
nlat = dims(1)
nlon = dims(2)
if (debug) then
;--- write to a netcdf file ---
fname = "output.nc"
system("/bin/rm -f "+fname)
fout = addfile(fname, "c")
;--- define mode on ---
setfileoption(fout, "DefineMode", True)
;--- define dimensions ---
dimNames = (/ "lat", "lon", "time" /)
dimSizes = (/ nlat, nlon, 1 /)
dimUnlim = (/ False, False, True /)
filedimdef(fout, dimNames, dimSizes, dimUnlim)
;--- define variables ---
filevardef(fout, "lat", typeof(lat1d), "lat")
filevardef(fout, "lon", typeof(lon1d), "lon")
filevardef(fout, "time", typeof(time), "time")
filevardef(fout, "uwnd", typeof(uwnd_rmap), (/ "time", "lat", "lon" /))
filevardef(fout, "vwnd", typeof(vwnd_rmap), (/ "time", "lat", "lon" /))
;--- add attributes ---
attr = True
attr at long_name = "Time"
attr at units = time at units
attr at calendar = time at calendar
filevarattdef(fout, "time", attr)
delete(attr)
attr = True
attr at long_name = "Longitude"
attr at units = "degrees_east"
filevarattdef(fout, "lon", attr)
delete(attr)
attr = True
attr at long_name = "Latitude"
attr at units = "degrees_north"
filevarattdef(fout, "lat", attr)
delete(attr)
;--- exit file definition mode ---
setfileoption(fout, "DefineMode", False)
;--- fill data ---
fout->lat = (/ lat1d /)
fout->lon = (/ lon1d /)
fout->time = (/ time /)
fout->uwnd = (/ uwnd_rmap /)
fout->vwnd = (/ vwnd_rmap /)
end if
;---------------------------------------------------------
; STEP (3)
; Write data in binary format
;---------------------------------------------------------
;--- create output file ---
fname = "wind.bin"
system("/bin/rm -f "+fname)
WAM::write_header(fname, llc(0), urc(0), llc(1), urc(1), 0.25, 0.25, nlat, nlon, 1)
do k = 0, ntime-1
print(k+" "+date_str(k))
u = uwnd_rmap(k,:,:)
v = vwnd_rmap(k,:,:)
WAM::write_uv(fname, date_str(k), nlon, nlat, u(lon|:,lat|:), v(lon|:,lat|:))
end do
end
More information about the ncl-talk
mailing list