[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