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