[gradsusr] Fortran to Read in GrADS Binary File
Jennifer M Adams
jadams21 at gmu.edu
Tue Jul 5 14:53:51 EDT 2016
The syntax of units field in a variable declaration depends on the data type: for grib is is a comma-delimited list of index numbers, for netcdf it is a comma-delimited list of coordinate axes, and for binary it is ignored.
—Jennifer
On 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>> 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>
_______________________________________________
gradsusr mailing list
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>
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> <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
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<mailto:gradsusr at gradsusr.org>
http://gradsusr.org/mailman/listinfo/gradsusr
--
Jennifer Miletta Adams
Center for Ocean-Land-Atmosphere Studies (COLA)
George Mason University
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20160705/da8bcd00/attachment-0001.html
More information about the gradsusr
mailing list