[gradsusr] fwrite gridding issue
Thomas Robinson
ter at hawaii.edu
Thu May 10 17:16:33 EDT 2012
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20120510/10eeefc9/attachment-0003.html
More information about the gradsusr
mailing list