[gradsusr] fwrite gridding issue
Eric Altshuler
ela at cola.iges.org
Thu May 10 20:20:23 EDT 2012
Tom,
I think I've found what's causing the problem. In your fortran code, you open unit 29 (the output from fwrite) as direct access, but the 'set fwrite' line in your grads script specifies sequential (-sq). This means that the control words in the fwrite output will be treated as data values when reading from unit 29, and the array 'chi' will not contain the correct data. To fix this, remove the -sq flag from your 'set fwrite' command.
I suppose it's unnecessary to write the output file (unit 291) as 48212 individual time steps, since your machine apparently has no problem with a 144*73*48212 array. I assume your DO loop on 70 is in response to my previous suggestion. If you wish, you can remove this loop and the declaration of chi_out, and revert to your original method of writing chi (which you've commented out).
In summary, your problem should go away if you do the following:
1. remove -sq flag from 'set fwrite'
2. just before fwrite, use 'set x 1 144' and 'set y 1 73' - leave your other 'set lon' and 'set lat' lines as they are
3. your ctl file should NOT have 'options sequential'
Eric
----- Original Message -----
From: "Thomas Robinson" <ter at hawaii.edu>
To: "GrADS Users Forum" <gradsusr at gradsusr.org>
Sent: Thursday, May 10, 2012 5:16:33 PM
Subject: Re: [gradsusr] fwrite gridding issue
I really appreciate all of our help. I'm just going to post all of my code now.
Here is my FORTRAN code as it stands now.
PROGRAM avg_chi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! This FORTRAN 77 program is used to calculate the daily mean values for the!!
!! entire data set. This mean value can then be used to calculate the !!
!! anomalies to be plotted. First, the data are opened in their individual !!
!! files. The full data set was divided into 30 files, so the parameter !!
!! files=32 represents those files. The numbers are converted to string by !!
!! writing their values to the string. !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!! Tom Robinson University of Hawaii 2012 ter at hawaii.edu !!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PARAMETER (nx=144) !grid points in x
PARAMETER (ny=73) !grid points in y
PARAMETER (tsy=1460) !Time Steps in a Year
PARAMETER (nt=709) !Number of time steps in each file
PARAMETER (nrec=709) !Number of records in each file
PARAMETER (files=68) !Number of input files
Parameter (rearth=6378100) !radius of the earth in meters
parameter (rpi=3.1415926) !pi
parameter(omega=7.292E-5) !Rotational speed of the earth
PARAMETER (south=-90) !Southern most latitude
PARAMETER (dy=111131)!change in y between a degree of lattitude at 45 deg
real chi (nx,ny,nt*files) ! Velocity potential
real chi_out(nx,ny)
real a_chi(nx,ny,tsy) !Average Velocity Potential
real clatrad (ny) !latitude in radians
real dx(iy) ! change in east-west direction (distance)
Integer x !Place holder for the file number
Character(2) xstring !Input file number > 9
Character(1) xstrLT1 !Input file number < 10
character(20) filename!Full input file name
it=1 !Initialize it
do 51 x=1,files!files !Loop through all of the file numbers
if (x.lt.10)then !Check to see if the number is below 10. If it
! is you have to use a single bit string
Write( xstrLT1, '(i1)' ) x !COnvert the integer to a string
write (6,*)"CHI_200hPa_",xstrLT1,'.dat'
filename="CHI_200hPa_"//xstrLT1//".dat" !Concatenate strings to
! get the file name
write (6,*)filename
!To open the file from a GrADS fwrite, you have to use direct
! access and specify the record length which is four times the
! dimensions
open (29, file=filename,
& status="OLD", access="direct",form="unformatted",
& recl=4*nx*ny)
do 52 irec=1,nrec !Go through each record of the input file
read(29,rec=irec)((chi(ix,iy,it),
& ix=1,nx),iy=1,ny)! Read in the data
it=it+1 !Update it
if (mod(it,100).eq.0)then
write(6,*)it
endif
52 continue
close (unit=29)
elseif (x.ge.10 .AND. x.lt.100)then !If the number of the file is
!GT 10, no 0 needs to be added
Write( xstring, '(i2)' ) x
write (6,*)"CHI_200hPa_",xstring,'.dat'
filename="CHI_200hPa_"//xstring//".dat" !Concatenate strings to
! get the file name
open (29, file=filename,
& status="OLD", access="direct",form="unformatted",
& recl=4*nx*ny)
do 53 irec=1,nrec
itime=it*x
read(29,rec=irec)((chi(ix,iy,it),
& ix=1,nx),iy=1,ny)
it=it+1
53 continue
close (unit=29)
endif
51 continue
!! Write out some numbers to see if they make sense
write (6,500)(chi(73,35,it)/1e6,it=5430,5436)
500 FORMAT(f7.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Open the output data file
open (unit=291,file="chi_33yr.dat", access="direct",
& status="unknown",recl=4*nx*ny)
!! Use chi_out to write each time step out as a record.
DO 70 itt=1,files*nt
do ix=1,nx
do iy=1,ny
chi_out(ix,iy)=chi(ix,iy,itt)
enddo
enddo
write (291, rec=itt) chi_out !Write each time step as a record
70 CONTINUE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Get the average values using the daily_ave subroutine
do it=1,tsy !Loop through all 1464 6 hour times in the year (leap year)
CALL daily_ave(chi,a_chi,it,nx,ny,nt,tsy)
if (mod(it,100).eq.0)then
write(6,*)it,a_chi(55,55,it)/1e6
endif
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Write the data out to a file
open (unit=290,file="mean_chi_33yr.dat", access="direct",
& status="unknown",recl=4*nx*ny*tsy)
write (290, rec=1) a_chi
! open (unit=291,file="chi_33yr.dat", access="direct",
! & status="unknown",recl=4*nx*ny*nt*files)
! write (291, rec=1) chi
stop
end
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
I left out the subroutines because I think they are working correctly and don't have anything to do with the problem I'm having since they are used below where the data is written out.
Here is the list of input files I have:
CHI_200hPa_10.dat CHI_200hPa_26.dat CHI_200hPa_41.dat CHI_200hPa_57.dat
CHI_200hPa_11.dat CHI_200hPa_27.dat CHI_200hPa_42.dat CHI_200hPa_58.dat
CHI_200hPa_12.dat CHI_200hPa_28.dat CHI_200hPa_43.dat CHI_200hPa_59.dat
CHI_200hPa_13.dat CHI_200hPa_29.dat CHI_200hPa_44.dat CHI_200hPa_5.dat
CHI_200hPa_14.dat CHI_200hPa_2.dat CHI_200hPa_45.dat CHI_200hPa_60.dat
CHI_200hPa_15.dat CHI_200hPa_30.dat CHI_200hPa_46.dat CHI_200hPa_61.dat
CHI_200hPa_16.dat CHI_200hPa_31.dat CHI_200hPa_47.dat CHI_200hPa_62.dat
CHI_200hPa_17.dat CHI_200hPa_32.dat CHI_200hPa_48.dat CHI_200hPa_63.dat
CHI_200hPa_18.dat CHI_200hPa_33.dat CHI_200hPa_49.dat CHI_200hPa_64.dat
CHI_200hPa_19.dat CHI_200hPa_34.dat CHI_200hPa_4.dat CHI_200hPa_65.dat
CHI_200hPa_1.dat CHI_200hPa_35.dat CHI_200hPa_50.dat CHI_200hPa_66.dat
CHI_200hPa_20.dat CHI_200hPa_36.dat CHI_200hPa_51.dat CHI_200hPa_67.dat
CHI_200hPa_21.dat CHI_200hPa_37.dat CHI_200hPa_52.dat CHI_200hPa_68.dat
CHI_200hPa_22.dat CHI_200hPa_38.dat CHI_200hPa_53.dat CHI_200hPa_6.dat
CHI_200hPa_23.dat CHI_200hPa_39.dat CHI_200hPa_54.dat CHI_200hPa_7.dat
CHI_200hPa_24.dat CHI_200hPa_3.dat CHI_200hPa_55.dat CHI_200hPa_8.dat
CHI_200hPa_25.dat CHI_200hPa_40.dat CHI_200hPa_56.dat CHI_200hPa_9.dat
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
Here is my GrADS script that I use to calculate the velocity potential and then write it out using fwrite:
function full_data(args)
*##############################################################################
*# Oh boy the fun stuff. This script opens a DODS file to get renalaysis ##
*# data that goes all the way back to 1979. It looks like it updates every ##
*# day (e'ryday e'ryday) so it should make thing interesting. The begin time##
*# is 00Z01JAN1979 and the time step is 360mn (6 hours). That's going to ##
*# complicate things. I suggest cutting it off at 18Z31DEC2011, that makes ##
*# the final time 48212. This is divisible by 68, 48212/68=709. We should ##
*# be able to do that many calculations (hopefully). We can! ##
*# ##
*# I took out the top 7 layers because the 200mb level is level 10. This ##
*# will make it run a little more efficient in theory. What I should do is ##
*# have it just do level 10. That would probably cut the time down a lot but##
*# I am a masocist. ##
*# This takes FOR-EV-ER! ##
*##############################################################################
*########### Tom Robinson University of Hawaii 2012 ter at hawaii.edu ############
*##############################################################################
*# Number of times divided - SHOULD BE 68 !!!
divnum=subwrd(args,1)
*# Number of file
filenum=subwrd(args,2)
say divnum ' ' filenum
*# Reinitialize everything, including the gxout
'reinit'
'set gxout contour'
'set looping off'
*# Open the internet file
'sdfopen http://nomad2.ncep.noaa.gov:9090/dods/reanalyses/reanalysis-2/6hr/pgb/pgb '
*#*****************************************************************************
*# SET UP GLOBAL VARIABLES AND MAP **
*#*****************************************************************************
*# get nt from the descriptor file
'q ctlinfo'
*tdef=sublin(result,9)
*nt=subwrd(tdef,2)
nt=48212
say nt
*# make sure looping is off. Who knows why it comes on sometimes :-(
'set looping off'
*# Set up the dimensions of the data. Values starting with n indicate the
*# number of data for the dimension that follows (x,y,z,t)
*dx=2.5
*dy=2.5
*#
nx=144
ny=73
nz=10
nt=nt
it=(((filenum-1)*nt)/divnum)+1
nt=filenum*nt/divnum
*# Set background to white
'set background 1'
*# Clear to get the white background
'clear'
'set grads off'
'set mpdset hires'
'set string 0'
*#* Set the map to black
'set map 0'
*#* Set the labels and title to black
'set annot 0'
*#radius on the earth in m
'define Re=6378100'
*# pi
'define pi=3.14159265'
*# ytom is the conversion of 1 degree to meters at 30N
'define ytom=110852.4248'
*# SET UP FOR ALL LEVELS AND TIMES **
*# set to full domain to define variables
'set t ' it ' ' nt
say 'set t ' it ' ' nt
'set lon 0 360'
'set lat -90 90'
*######## Z is set up for 200mb ONLY!
'set lev 200'
*#*****************************************************************************
*# CALCULATE VELOCITY POTENTIAL (chi) **
*#*****************************************************************************
*#** VELOCITY POTENTIAL FOR OPENGRADS ONLY ****
*# Calculate the velocity potential using the opengrads function fish_chi
'define chi = fish_chi(UGRDprs,VGRDprs)'
*******************************************************************************
** DISPLAY CHI AND DIV **
*#*****************************************************************************
*# Set up colors 16-20 are green, 21-25 are blue
gnum = 16
bnum = 21
count = 1
cnum = 25
while (count<=5)
*# Put the rgb colors in using gnum as gree and bnum as blue
'set rgb 'gnum' 0 ' cnum ' 0'
'set rgb 'bnum' 0 0 ' cnum
*# Increment the counter and color numbers
count=count+1
gnum=gnum+1
bnum=bnum+1
*# Increase the rgb number by 55
cnum=cnum+55
endwhile
'set gxout shaded'
*# Set up the levels and colors using blue for negative and green for positive
'set clevs -20 -15 -10 -5 0 5 10 15 20'
'set ccols 22 23 24 25 1 1 20 19 18 17 '
*# Change to a tropical latitude domain for display purposes
'set lat -30 30'
'set lon 0 360'
'set t 'it
'set lev 200 200'
*# Display the velocity potential and divide by a million
'd (chi/1e6)'
*# Display the color bar
'cbarn_black'
*# Draw the title with the date
'q time'
date=subwrd(result,3)
'q dim'
rec=sublin(result,4)
level=subwrd(rec,6)
'draw title 'level'hPa Velocity Potential 'date
** get the month
month=substr(date,6,3)
** get the year
year=substr(date,9,4)
** get the day
day=substr(date,4,2)
** get the hour
hour=substr(date,1,3)
say month ' ' year
'printim 'year''month''day''hour'_'level'hPa.jpg'
'printim 'year'_'filenum'.jpg'
*#*****************************************************************************
*# WRITE VOLOCITY POTENTIAL OUT TO A BINARY FILE **
*#*****************************************************************************
*# set up global domain once again
'set looping off'
'set lev 200 '
'set lat -90 90'
'set lon 0 360'
*#**************************************************************
say 'write data out to binary file called CHI.200hPa'filenum'.dat'
'set fwrite -le -sq -cl CHI_200hPa_'filenum'.dat'
'set gxout fwrite'
while (it <= nt)
'set t 'it
'display chi'
it=it+1
endwhile
'disable fwrite'
*******************************************************************************
*set fwrite <-be or -le> <-sq or -st> <-ap or -cl> fname
*Sets the filename for data output as well as byte ordering and data format.
* fname output filename (default = grads.fwrite)
* -be output data byte ordering is big endian
* -le output data byte ordering is little endian
* -sq output data format is sequential
* -st output data format is stream (default)
* -ap output data is appended to existing file
* -cl output data replaces existing file if it exists (default)
*******************************************************************************
'set looping off'
** Looping should be turned off
'set gxout contour'
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
On Thu, May 10, 2012 at 10:56 AM, Eric Altshuler < ela at cola.iges.org > wrote:
Tom,
There are three things I'll point out here:
1. You seem to be trying to write out all the data at once, in a huge array dimensioned (144,73,48212). This array may require more memory than your computer has. Instead, try declaring your array as (144,73) and set up a DO loop in your fortran program to write out each time step individually. You'll need to substantially modify how your program reads the input data, too. If you send your complete fortran code, I might be able to modify it to implement my suggested method (no guarantees, though). This will eliminate the need to break up the 48212 time steps into 68*709.
2. What is the size of the input files to your fortran program (i.e. the files written by fwrite)? This information will indicate if you have a wraparound problem.
3. Your ctl file has the line:
chi 0 x,y,t 1 ** velocity potential
Since your data is binary, I think this line should be:
chi 0 99 ** velocity potential
Eric
----- Original Message -----
From: "Thomas Robinson" < ter at hawaii.edu >
To: "GrADS Users Forum" < gradsusr at gradsusr.org >
Sent: Thursday, May 10, 2012 2:03:09 PM
Subject: Re: [gradsusr] fwrite gridding issue
The file is the correct size. 144*73*48212*4= 2027218176
-rw-r--r-- 1 ter businger 2027218176 May 9 15:57 chi_33yr.dat
My computer crashes if I try to do all 48212, so I had to break it up in 68 separate files each with 709 time steps. So:
PARAMETER (nt=709) !Number of time steps in each file
PARAMETER (nrec=709) !Number of records in each file
PARAMETER (files=68) !Number of input files
The FORTRAN program reads in all of the files and stores the data in one variable chi, then it goes written out to the file chi_33yr.dat.
My FORTRAN program also calculates a 33 year 6-hourly mean, but the variable chi isn't affected by that and that is written to a separate file which is also having the same problem.
-Tom
On Thu, May 10, 2012 at 7:51 AM, Jennifer Adams < jma at cola.iges.org > wrote:
I'm no fortran expert, but I think "direct" means the 4-byte headers and footers for each record are not written, so you should remove 'options sequential'. According to your descriptor, your data file (chi_33yr.dat) should be 144*73*48212*4 bytes large. Is it? What is 'files'?
--Jennifer
On May 10, 2012, at 1:44 PM, Thomas Robinson wrote:
Here is my code that opens and writes out the data:
open (unit=291,file="chi_33yr.dat", access="direct",
& status="unknown",recl=4*nx*ny*nt*files)
write (291, rec=1) chi
The variable chi has dimensions nx,ny,nt*files .
I added options sequential to the descriptor file and the same thing is happening. Here is my descriptor file:
dset ^chi_33yr.dat
title Velocity Potential
options sequential
undef 9.999e+20
xdef 144 linear 0 2.5
ydef 73 linear -90 2.5
zdef 1 levels 200
tdef 48212 linear 00Z01JAN1979 360mn
vars 1
chi 0 x,y,t 1 ** velocity potential
endvars
All of your input is very much appreciated.
-Tom
On Fri, May 11, 2012 at 3:16 AM, Jennifer Adams < jma at cola.iges.org > wrote:
Tom, Did you include 'options sequential' in your descriptor file for your new binary file? --Jennifer
On May 10, 2012, at 8:49 AM, Raghu Reddy wrote:
I am not a grads user, just responding to the FORTRAN/C issue.
If you do a normal sequential FORTRAN write, it writes an end-of-record marker at the beginning * and * at the end of * every * write. Since you mentioned that the data is skewed, I assume you wrote all the 48,000 values with one write statement, and that would result in one EOR marker at the beginning and one EOR marker at the end.
If you don’t want these EOR markers you could write the file out as direct access files in Fortran.
Hope this helps.
Thanks,
--Raghu
From: gradsusr-bounces at gradsusr.org [mailto: gradsusr-bounces at gradsusr.org ] On Behalf Of Thomas Robinson
Sent: Wednesday, May 09, 2012 11:36 PM
To: GrADS Users Forum
Subject: [gradsusr] fwrite gridding issue
Aloha,
I used fwrite to write out to a binary file in order to create a large file with about 48000 time steps for global data (because I am crazy). I noticed that when I use fortran to create a new data file that has all 48000, it doesn't match with what the original data plotted, instead it seems skewed. No manipulations were done to the data, it was just read it all in and write it out. When I tried to loop the data, I further noticed that the northern border looped around the south and then moved upwards back to the north. It just keeps looping around and it's odd because the top boundary of the map isn't continuous with the bottom boundary (unlike the side boundaries which are continuous). Anyways, I am wondering why this happened and what I can do to make it work out correctly.
Here is my fwrite code. it is the first time step, nt is the last time step, chi is the velocity potential that was calculated using 'define chi = fish_chi(UGRDprs,VGRDprs)' :
*# set up global domain
'set lev 200 '
'set lat -90 90'
'set lon 0 360'
*#**************************************************************
say 'write data out to binary file called CHI.200hPa'filenum'.dat'
'set fwrite -le -sq -cl CHI_200hPa_'filenum'.dat'
'set gxout fwrite'
while (it <= nt)
'set t 'it
'display chi'
it=it+1
endwhile
'disable fwrite'
Mahalo for your help!
-Tom
--
Tom Robinson
President Graduate Student Organization
Student Caucus Representative for the Graduate Student Organization
Graduate Student - Department of Meteorology
732-718-2323 _______________________________________________
gradsusr mailing list
gradsusr at gradsusr.org
http://gradsusr.org/mailman/listinfo/gradsusr
--
Jennifer M. Adams
IGES/COLA
4041 Powder Mill Road, Suite 302
Calverton, MD 20705
jma at cola.iges.org
_______________________________________________
gradsusr mailing list
gradsusr at gradsusr.org
http://gradsusr.org/mailman/listinfo/gradsusr
--
Tom Robinson
President Graduate Student Organization
Student Caucus Representative for the Graduate Student Organization
Graduate Student - Department of Meteorology
732-718-2323
_______________________________________________
gradsusr mailing list
gradsusr at gradsusr.org
http://gradsusr.org/mailman/listinfo/gradsusr
--
Jennifer M. Adams
IGES/COLA
4041 Powder Mill Road, Suite 302
Calverton, MD 20705
jma at cola.iges.org
_______________________________________________
gradsusr mailing list
gradsusr at gradsusr.org
http://gradsusr.org/mailman/listinfo/gradsusr
--
Tom Robinson
President Graduate Student Organization
Student Caucus Representative for the Graduate Student Organization
Graduate Student - Department of Meteorology
732-718-2323
_______________________________________________
gradsusr mailing list
gradsusr at gradsusr.org
http://gradsusr.org/mailman/listinfo/gradsusr
_______________________________________________
gradsusr mailing list
gradsusr at gradsusr.org
http://gradsusr.org/mailman/listinfo/gradsusr
--
Tom Robinson
President Graduate Student Organization
Student Caucus Representative for the Graduate Student Organization
Graduate Student - Department of Meteorology
732-718-2323
_______________________________________________
gradsusr mailing list
gradsusr at gradsusr.org
http://gradsusr.org/mailman/listinfo/gradsusr
More information about the gradsusr
mailing list