[gradsusr] Fortran to Read in GrADS Binary File
Lydia Rill
lydia.d.rill at gmail.com
Wed Jul 6 09:42:25 EDT 2016
Chuck- I printed out the data like you suggested, and almost all the values
from the tmax(i,j) array are missing while the tmax_all(i,j) array has
realistic values. I decided to read the whole line in first and then write
it out by x,y and it seems to be working!
Jim- I double checked the file sizes and they are correct. Thanks for that
useful tip!
Here is the final script I used and it is giving me realistic data in the
correct grid box!
integer irec,i,j,t,days
character(len=100) infile,outfile
real tmax_all(464,224), tmin_all(464,224),prcp_all(464,224),srad_all(464,224
)
infile='/Volumes/WebData/NLDAS2/DailyNLDAS/daily.bin'
outfile='/Volumes/WebData/NLDAS2/Fortran/testdaily.bin'
open(17,file=infile,form="unformatted",access="direct",recl=464*224*4,
status="old")
open(18,file=outfile,form="unformatted",access="direct",recl=32*32*4,status=
"unknown")
days=10
irec=1
do t=1,days
read(17,REC=irec) tmax_all
read(17,REC=irec+1) tmin_all
read(17,REC=irec+2) prcp_all
read(17,REC=irec+3) srad_all
write(18,REC=irec) ((tmax_all(i,j),i=129,160),j=129,160)
write(18,REC=irec+1) ((tmin_all(i,j),i=129,160),j=129,160)
write(18,REC=irec+2) ((prcp_all(i,j),i=129,160),j=129,160)
write(18,REC=irec+3) ((srad_all(i,j),i=129,160),j=129,160)
irec=irec+4
enddo
Thank you all for your help! I really appreciate it.
Lydia
On Tue, Jul 5, 2016 at 3:57 PM, Charles Seman - NOAA Federal <
charles.seman at noaa.gov> wrote:
> Lydia,
>
> It may be better to put the tmax_all and tmax printout on the same line:
>
> print *, 't, i, j, tmax_all(i,j), tmax(i,j) : ',t, i, j, tmax_all(i,j),
> tmax(i,j)
>
> this should make it easier to compare the values...
>
> you could also try printing the difference: tmax_all(i,j) - tmax(i,j)
>
> Chuck
>
> On 07/05/2016 03:36 PM, Charles Seman - NOAA Federal wrote:
> > Lydia,
> >
> > I would print directly in your Fortran program, within the "t" do-loop.
> > Try something like:
> > ...
> > real tmax_all(464,224)
> > ...
> > do t=1,days
> > read(17,REC=irec) ((tmax(i,j),i=129,160),j=129,160)
> > read(17,REC=irec) tmax_all ! to keep things simple, read in the
> > entire array
> > i=129
> > do j=129,160
> > print *, 't, i, j, tmax_all(i,j) : ',t, i, j, tmax_all(i,j)
> > print *, 't, i, j, tmax(i,j) : ',t, i, j, tmax(i,j)
> > enddo
> > j=129
> > do i=129,160
> > print *, 't, i, j, tmax_all(i,j) : ',t, i, j, tmax_all(i,j)
> > print *, 't, i, j, tmax(i,j) : ',t, i, j, tmax(i,j)
> > enddo
> > i=160
> > do j=129,160
> > print *, 't, i, j, tmax_all(i,j) : ',t, i, j, tmax_all(i,j)
> > print *, 't, i, j, tmax(i,j) : ',t, i, j, tmax(i,j)
> > enddo
> > j=160
> > do i=129,160
> > print *, 't, i, j, tmax_all(i,j) : ',t, i, j, tmax_all(i,j)
> > print *, 't, i, j, tmax(i,j) : ',t, i, j, tmax(i,j)
> > enddo
> > stop
> > ...
> > enddo
> >
> > The above printouts are designed to print the 4 edges of the desired
> > 32x32 array of data, and may reveal something. You could also try
> > printing across the middle of the array (horizontally and vertically).
> >
> > If nothing is amiss, then try reading and writing just one variable to a
> > GrADS binary file, and adjust ctl file... see if this looks OK... if so,
> > then try another variable...
> >
> > Hope this helps,
> > Chuck
> >
> > On 07/05/2016 02:56 PM, Lydia Rill wrote:
> >> Hi Chuck,
> >>
> >> This is my first time using GrADS so I'm still pretty fuzzy on the
> >> specifics in the control file. My understanding is the 1st number is the
> >> vertical level, and the 2nd value is an ID label. The 3rd value I am not
> >> sure, and the 4th value is the height at which the variable was
> >> measured. I took these values directly from a control file created using
> >> the grib2ctl function when processing the raw data which were grib
> >> files. (I processed them by aggregating hourly to daily data).
> >>
> >> When you say print out a row and column, am I able to print the data as
> >> text in Fortran? Or should I use GrADS?
> >>
> >> Thanks,
> >>
> >> Lydia
> >>
> >>
> >>
> >> On Tue, Jul 5, 2016 at 2:31 PM, Charles Seman - NOAA Federal
> >> <charles.seman at noaa.gov <mailto:charles.seman at noaa.gov>> wrote:
> >>
> >> Lydia,
> >>
> >> Another thing: I don't understand the characters like "12,105,2
> >> **" in
> >> the vars records in your ctl file:
> >>
> >> TMAX 0 12,105,2 ** maximum temperature at 2m above surface [C]
> >>
> >> TMIN 0 11,105,2 ** minimum temperature at 2m above surface [C]
> >>
> >> PRCP 0 61,1,0 ** total precipitation backward accumulated [kg/m^2
> >> or mm]
> >>
> >> SRAD 0 204,1,0 ** total surface downward shortwave radiation flux
> >> [mJ/m^2]
> >>
> >> How are these characters used by GrADS? I would expect something
> >> like
> >> "0,y,x"
> >>
> >> Chuck
> >>
> >> On 07/05/2016 02:21 PM, Charles Seman - NOAA Federal wrote:
> >> > Lydia,
> >> >
> >> > I would try to verify some of the values in your arrays... To
> >> start,
> >> > define a tmax_all(464,224), and read one record into it:
> >> > read(17,REC=irec) tmax_all. Then, print out a common row and a
> >> column
> >> > from it and from tmax to verify you get the same values. Do
> >> this for
> >> > just one time level, and stop. This is slow and tedious work.
> >> Unless
> >> > you have access to TotalView or some other visual debugger, then
> >> you can
> >> > step through your code line by line and should be able to find it.
> >> >
> >> > Hope this helps,
> >> > Chuck
> >> >
> >> > On 07/05/2016 01:41 PM, Lydia Rill wrote:
> >> >> Thanks Chuck! I double checked the record length and it is
> >> indeed 4.
> >> >>
> >> >> Thank you Jim! Sorry my code was a little messy, I have tried
> many
> >> >> different methods.
> >> >> Good point about the lower left hand corner, that was one of the
> >> issues
> >> >> here.
> >> >>
> >> >> When I use 32 x steps, 32 y steps, and 10 time steps I use a
> >> record
> >> >> length of 4096 for the output file and I get a file size of 160K,
> >> which
> >> >> roughly matches 32*32*10*4*4 = 16385
> >> >>
> >> >> I am now getting realistic numbers for my output! However, when I
> >> >> display the data in GrADS for one time I get only parts of the
> 4x4
> >> >> degree box, by latitude, for each variable, instead of getting
> >> the data
> >> >> for each variable for the whole box. For example in the middle
> >> of the
> >> >> box I get tmax for 43.5-44N, tmin for 43-43.5N, prcp around
> >> 43N, and
> >> >> srad around 42.5N. Any ideas on why this is happening?
> >> >>
> >> >> Here is the updated code:
> >> >>
> >> >> integer irec,i,j,t,days
> >> >>
> >> >> character(len=100) infile,outfile
> >> >>
> >> >> real tmin(464,224),tmax(464,224),prcp(464,224),srad(464,224)
> >> >>
> >> >> infile='/Volumes/WebData/NLDAS2/DailyNLDAS/daily.bin'
> >> >>
> >> >> outfile='/Volumes/WebData/NLDAS2/Fortran/testdaily.bin'
> >> >>
> >> >>
> >>
> >>
> open(17,file=infile,form="unformatted",access="direct",recl=32*32*4,status="old")
> >>
> >> >>
> >> >>
> >> >>
> >>
> >>
> open(18,file=outfile,form="unformatted",access="direct",recl=32*32*4,status="unknown")
> >>
> >> >>
> >> >>
> >> >>
> >> >> days=10
> >> >>
> >> >> irec=1
> >> >>
> >> >> do t=1,days
> >> >>
> >> >> read(17,REC=irec) ((tmax(i,j),i=129,160),j=129,160)
> >> >>
> >> >> read(17,REC=irec+1) ((tmin(i,j),i=129,160),j=129,160)
> >> >>
> >> >> read(17,REC=irec+2) ((prcp(i,j),i=129,160),j=129,160)
> >> >>
> >> >> read(17,REC=irec+3) ((srad(i,j),i=129,160),j=129,160)
> >> >>
> >> >> write(18,REC=irec) ((tmax(i,j),i=129,160),j=129,160)
> >> >>
> >> >> write(18,REC=irec+1) ((tmin(i,j),i=129,160),j=129,160)
> >> >>
> >> >> write(18,REC=irec+2) ((prcp(i,j),i=129,160),j=129,160)
> >> >>
> >> >> write(18,REC=irec+3) ((srad(i,j),i=129,160),j=129,160)
> >> >>
> >> >> irec=irec+4
> >> >>
> >> >> enddo
> >> >>
> >> >>
> >> >> Here is the updated control file:
> >> >>
> >> >> dset ^testdaily.bin
> >> >>
> >> >> undef 9.999E+20
> >> >>
> >> >> title Daily NLDAS2 Data for 04Z through 03Z
> >> >>
> >> >> xdef 32 linear -108.9375 0.125
> >> >>
> >> >> ydef 32 linear 41.0625 0.125
> >> >>
> >> >> tdef 10 linear 00Z02Jan1979 1dy
> >> >>
> >> >> zdef 1 linear 1 1
> >> >>
> >> >> vars 4
> >> >>
> >> >> TMAX 0 12,105,2 ** maximum temperature at 2m above surface [C]
> >> >>
> >> >> TMIN 0 11,105,2 ** minimum temperature at 2m above surface [C]
> >> >>
> >> >> PRCP 0 61,1,0 ** total precipitation backward accumulated [kg/m^2
> >> or mm]
> >> >>
> >> >> SRAD 0 204,1,0 ** total surface downward shortwave radiation flux
> >> >> [mJ/m^2]
> >> >>
> >> >> ENDVARS
> >> >>
> >> >>
> >> >> Thank you both so much for your help,
> >> >>
> >> >> Lydia
> >> >>
> >> >> On Tue, Jul 5, 2016 at 12:44 PM, James T. Potemra
> >> <jimp at hawaii.edu <mailto:jimp at hawaii.edu>
> >> >> <mailto:jimp at hawaii.edu <mailto:jimp at hawaii.edu>>> wrote:
> >> >>
> >> >> Hi Lydia,
> >> >>
> >> >> You have a couple problems. First, you are reading in tmax,
> >> tmin,
> >> >> prcp then srad, but then you write out tmin, tmin prcp and
> >> srad
> >> >> (tmin twice). Second, your program is reading and writing
> >> to the
> >> >> same variables, and then only the first 8x8 (lower left
> >> "corner" of
> >> >> the data sets). Perhaps these are all missing ?
> >> >>
> >> >> Finally, the record length is defined to be 464 by 224,
> >> then you
> >> >> only define (via read) and write the first 8 by 8, which
> >> should work
> >> >> ok, but recognize that the rest of the 464x224 array has
> >> not been
> >> >> defined. So, if you want to use your control file with "xdef
> >> 8" and
> >> >> "ydef 8" you'll need to change your "open" statement for the
> >> output
> >> >> file to be 8*8*4 and not 464*224*4. Alternately you could
> >> change
> >> >> the xdef/ydef lines.
> >> >>
> >> >> Note as a quick check, your output file should be 8*8*10*4*4
> >> >> (x*y*t*4) bytes large for that control file.
> >> >>
> >> >> Jim
> >> >>
> >> >> On 7/5/16 3:57 AM, Lydia Rill wrote:
> >> >>> Hi all,
> >> >>>
> >> >>> I am a new grads user and I have never worked with binary
> >> files
> >> >>> before. I have a large binary file I created in grads (code
> >> shown
> >> >>> below). I want to split up the binary file by latitude and
> >> >>> longitude degree into smaller binary files. While this is
> >> possible
> >> >>> in GrADS, it is too slow for my needs, so I am trying to do
> >> it in
> >> >>> Fortran 90.
> >> >>>
> >> >>> However, I am having problems reading the GrADS binary
> >> file in
> >> >>> Fortran. It gives me the missing value for each of my
> >> variables at
> >> >>> every time step and every location.
> >> >>>
> >> >>> The binary file contains NLDAS 2 forcing data, with 464 x
> >> steps,
> >> >>> 224 y steps, and 13688 daily time steps for 4 variables
> >> (only 1 Z
> >> >>> level). It is a large file at 22GB, and I can successfully
> >> display
> >> >>> the variables in GrADS from this binary file. This binary
> >> file was
> >> >>> created from data in Grib files.
> >> >>>
> >> >>> I have tried many various ways of reading in the binary
> data,
> >> >>> using x and y loops, or time loops, or looping through the 4
> >> >>> variables, but nothing has been successful.
> >> >>>
> >> >>> I am working on a Mac (El Capitan).
> >> >>>
> >> >>> Here is the GrADS script I used to create the binary file:
> >> >>>
> >> >>> functionmain(args)
> >> >>>
> >> >>> counthr=subwrd(args,1)
> >> >>>
> >> >>> 'reinit'
> >> >>>
> >> >>> 'open /Volumes/WebData/NLDAS2/ProcessingFiles/24hours.ctl'
> >> >>>
> >> >>> 'set undef 9.999e20'
> >> >>>
> >> >>> 'set gxout fwrite'
> >> >>>
> >> >>> 'set fwrite -ap
> /Volumes/WebData/NLDAS2/DailyNLDAS/daily.bin'
> >> >>>
> >> >>>
> >> >>> 'set t 1'
> >> >>>
> >> >>> t1=1+4
> >> >>>
> >> >>> t1p=1+5
> >> >>>
> >> >>> while(t1<=counthr+5-23)
> >> >>>
> >> >>> t2=t1+23
> >> >>>
> >> >>> t2p=t1p+23
> >> >>>
> >> >>> 'd (max(TMP2m,t='t1',t='t2')-273.15)'
> >> >>>
> >> >>> 'd (min(TMP2m,t='t1',t='t2')-273.15)'
> >> >>>
> >> >>> 'd sum(APCPsfc,t='t1p',t='t2p')'
> >> >>>
> >> >>> 'd (ave(DSWRFsfc,t='t1',t='t2')*24*3600/1000000)'
> >> >>>
> >> >>> t1=t1+24
> >> >>>
> >> >>> t1p=t1p+24
> >> >>>
> >> >>> endwhile
> >> >>>
> >> >>> 'quit'
> >> >>>
> >> >>>
> >> >>> Here is the Fortran script to read the binary file
> >> >>>
> >> >>> ProgramReadbin
> >> >>>
> >> >>> implicitnone
> >> >>>
> >> >>> integer irec,i,j,t,days
> >> >>>
> >> >>> real tmin(464,224),tmax(464,224),prcp(464,224),srad(464,224)
> >> >>>
> >> >>>
> OPEN(17,file='/Volumes/WebData/NLDAS2/DailyNLDAS/daily.bin',&
> >> >>>
> >> >>>
> >> &STATUS='Old',FORM='UNFORMATTED',ACCESS='DIRECT',RECL=464*224*4)
> >> >>>
> >> >>>
> >> >>>
> >> OPEN(18,file='/Volumes/WebData/NLDAS2/DailyNLDAS/testdaily.bin',&
> >> >>>
> >> >>>
> >>
> &STATUS='UNKNOWN',FORM='UNFORMATTED',ACCESS='DIRECT',RECL=464*224*4)
> >> >>>
> >> >>> days=10
> >> >>>
> >> >>> irec=1
> >> >>>
> >> >>> do t=1,days
> >> >>>
> >> >>> read(17,REC=irec) ((tmax(i,j),i=1,8),j=1,8)
> >> >>>
> >> >>> read(17,REC=irec+1) ((tmin(i,j),i=1,8),j=1,8)
> >> >>>
> >> >>> read(17,REC=irec+2) ((prcp(i,j),i=1,8),j=1,8)
> >> >>>
> >> >>> read(17,REC=irec+3) ((srad(i,j),i=1,8),j=1,8)
> >> >>>
> >> >>> write(18,REC=irec) ((tmin(i,j),i=1,8),j=1,8)
> >> >>>
> >> >>> write(18,REC=irec+1) ((tmin(i,j),i=1,8),j=1,8)
> >> >>>
> >> >>> write(18,REC=irec+2) ((prcp(i,j),i=1,8),j=1,8)
> >> >>>
> >> >>> write(18,REC=irec+3) ((srad(i,j),i=1,8),j=1,8)
> >> >>>
> >> >>> irec=irec+4
> >> >>>
> >> >>> enddo
> >> >>>
> >> >>> end program
> >> >>>
> >> >>>
> >> >>> The corresponding control file is:
> >> >>>
> >> >>> set ^testdaily.bin
> >> >>>
> >> >>> undef 9.999E+20
> >> >>>
> >> >>> title Daily NLDAS2 Data for00Z through 23Z
> >> >>>
> >> >>> xdef 8linear -124.93750.125
> >> >>>
> >> >>> ydef 8linear 25.06250.125
> >> >>>
> >> >>> tdef 10linear 00Z02Jan1979 1dy
> >> >>>
> >> >>> zdef 1linear 11
> >> >>>
> >> >>> vars 4
> >> >>>
> >> >>> TMAX 012,105,2**maximum temperature at 2m above surface [C]
> >> >>>
> >> >>> TMIN 011,105,2**minimum temperature at 2m above surface [C]
> >> >>>
> >> >>> PRCP 061,1,0**total precipitation backward accumulated
> >> [kg/m^2ormm]
> >> >>>
> >> >>> SRAD 0204,1,0**total surface downward shortwave radiation
> >> flux
> >> >>> [mJ/m^2]
> >> >>>
> >> >>> ENDVARS
> >> >>>
> >> >>>
> >> >>>
> >> >>> Thanks,
> >> >>>
> >> >>> Lydia
> >> >>>
> >> >>> --
> >> >>> Lydia D. Rill
> >> >>> M.S., Michigan State University 2016
> >> >>> B.S., Valparaiso University 2014
> >> >>> (224) 406-5130
> >> >>> Lydia.D.Rill at gmail.com <mailto:Lydia.D.Rill at gmail.com>
> >> <mailto:Lydia.D.Rill at gmail.com <mailto:Lydia.D.Rill at gmail.com>>
> >> >>>
> >> >>>
> >> >>> _______________________________________________
> >> >>> gradsusr mailing list
> >> >>> gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
> >> <mailto:gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>>
> >> >>> http://gradsusr.org/mailman/listinfo/gradsusr
> >> >>
> >> >>
> >> >> _______________________________________________
> >> >> gradsusr mailing list
> >> >> gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
> >> <mailto:gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>>
> >> >> http://gradsusr.org/mailman/listinfo/gradsusr
> >> >>
> >> >>
> >> >>
> >> >>
> >> >> --
> >> >> Lydia D. Rill
> >> >> M.S., Michigan State University 2016
> >> >> B.S., Valparaiso University 2014
> >> >> (224) 406-5130 <tel:%28224%29%20406-5130>
> >> >> Lydia.D.Rill at gmail.com <mailto:Lydia.D.Rill at gmail.com>
> >> <mailto:Lydia.D.Rill at gmail.com <mailto:Lydia.D.Rill at gmail.com>>
> >> >>
> >> >>
> >> >> _______________________________________________
> >> >> gradsusr mailing list
> >> >> gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
> >> >> http://gradsusr.org/mailman/listinfo/gradsusr
> >> >>
> >> >
> >>
> >> --
> >>
> >> Please note that Charles.Seman at noaa.gov
> >> <mailto:Charles.Seman at noaa.gov> should be considered my NOAA
> >> email address, not cjs at gfdl.noaa.gov <mailto:cjs at gfdl.noaa.gov>.
> >>
> >> ********************************************************************
> >> Charles Seman
> >> Charles.Seman at noaa.gov <mailto:Charles.Seman at noaa.gov>
> >> U.S. Department of Commerce / NOAA / OAR
> >> Geophysical Fluid Dynamics Laboratory voice: (609)
> >> 452-6547 <tel:%28609%29%20452-6547>
> >> 201 Forrestal Road fax: (609)
> >> 987-5063 <tel:%28609%29%20987-5063>
> >> Princeton, NJ 08540-6649
> http://www.gfdl.noaa.gov/~cjs/
> >> ********************************************************************
> >>
> >> "The contents of this message are mine personally and do not
> >> reflect any
> >> official or unofficial position of the United States Federal
> >> Government,
> >> the United States Department of Commerce, or NOAA."
> >> _______________________________________________
> >> gradsusr mailing list
> >> gradsusr at gradsusr.org <mailto:gradsusr at gradsusr.org>
> >> http://gradsusr.org/mailman/listinfo/gradsusr
> >>
> >>
> >>
> >>
> >> --
> >> Lydia D. Rill
> >> M.S., Michigan State University 2016
> >> B.S., Valparaiso University 2014
> >> (224) 406-5130
> >> Lydia.D.Rill at gmail.com <mailto:Lydia.D.Rill at gmail.com>
> >>
> >>
> >> _______________________________________________
> >> gradsusr mailing list
> >> gradsusr at gradsusr.org
> >> http://gradsusr.org/mailman/listinfo/gradsusr
> >>
> >
>
> --
>
> Please note that Charles.Seman at noaa.gov should be considered my NOAA
> email address, not cjs at gfdl.noaa.gov.
>
> ********************************************************************
> Charles Seman Charles.Seman at noaa.gov
> U.S. Department of Commerce / NOAA / OAR
> Geophysical Fluid Dynamics Laboratory voice: (609) 452-6547
> 201 Forrestal Road fax: (609) 987-5063
> Princeton, NJ 08540-6649 http://www.gfdl.noaa.gov/~cjs/
> ********************************************************************
>
> "The contents of this message are mine personally and do not reflect any
> official or unofficial position of the United States Federal Government,
> the United States Department of Commerce, or NOAA."
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
>
--
Lydia D. Rill
M.S., Michigan State University 2016
B.S., Valparaiso University 2014
(224) 406-5130
Lydia.D.Rill at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20160706/0506d424/attachment-0001.html
More information about the gradsusr
mailing list