[gradsusr] Fortran to Read in GrADS Binary File
Charles Seman - NOAA Federal
charles.seman at noaa.gov
Tue Jul 5 14:21:19 EDT 2016
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>> 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>
>>
>>
>> _______________________________________________
>> gradsusr mailing list
>> 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>
> 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."
More information about the gradsusr
mailing list