Multi-level Station Data
Brian Reen
snincr at YAHOO.COM
Fri Sep 16 09:16:20 EDT 2005
GrADS users,
I've determined the errors in my previous attempt to
use multi-level station data.
1) There should be a single write statement for all of
the data for a given station at a given time.
2) In checking that the data has been correctly
written using the method I previously indicated ('set
gxout stat', 'set stnprint on', 'd a'), one must set
the vertical coordinate to one where data exists (e.g.
'set z 15').
Therefore, the following Fortran code and .ctl file
are an example of very simple multi-level station
data. To use this code do the following:
1) Compile and run the Fortran code.
2) stnmap -i testcase.ctl
3) Open GrADS and issue the following commands:
a) open testcase.ctl
b) set gxout stat
c) set stnprint on
d) set z 15
e) d a
This should result in output that ends in:
Printing station values: #obs = 1
OB ID LON LAT LEV VAL
1 35 45 15 12
***************************************************
Begin Fortran Code
***************************************************
PROGRAM TESTCASE
IMPLICIT NONE
integer :: outunit = 12 !Output file unit
character(len=8) :: stationid !Station identifier
for grads file
real :: lat,lon !Latitude, longitude
real :: timeshift !Time shift of this ob
integer :: nlev !Number of vertical layers (+1 if
flag==1)
integer :: flag !Is there sfc-only data
real :: height(2) !Height of ob
real :: ob(2) !Observed value
integer :: levelcount !Counter over levels
!Open the output file
open(outunit,file='testcase.dat',status='replace',action='write',&
form='unformatted',access='sequential')
write(stationid,'I8'),1
lat = 45.0
lon = 35.0
timeshift = 0.0
nlev = 2
flag = 0
!Write header
write(outunit),stationid,lat,lon,timeshift,nlev,flag
!Write data
height(1) = 15.0
height(2) = 25.0
ob(1) = 12.0
ob(2) = 22.0
write(outunit),(height(levelcount),ob(levelcount),levelcount=1,2)
!Write time terminator
write(outunit),'TIMETERM',0.0,0.0,0.0,0,0
close(outunit)
END PROGRAM TESTCASE
***************************************************
End Fortran Code
***************************************************
***************************************************
Begin .ctl file
***************************************************
DSET ^testcase.dat
DTYPE station
STNMAP testcase.map
UNDEF -999.0
TITLE Test Data
OPTIONS SEQUENTIAL
TDEF 1 linear 00:00Z01may2002 01mn
VARS 1
a 1 99 no description
ENDVARS
***************************************************
End .ctl file
***************************************************
To illustrate how to write multi-level station data, I
also have code that attempts to create multi-level
station data with an arbitrary number of surface
variables, multi-level variables, stations, times, and
levels. That code follows:
PROGRAM TSMIN
IMPLICIT NONE
real, allocatable, dimension(:) :: lat,lon
!Latitude/longitude of stations
character(len=8), allocatable, dimension(:) ::
stationid !Station identifier for grads file
integer :: outunit = 12 !Output file unit
integer :: timecounter !Counter over output times
integer :: timecountermax !Total number of times
integer :: stationcounter !Counter over stations
integer :: stationcountermax !Total number of
stations
integer :: levelcounter !Counter over vertical
levels
integer :: levelcountermax !Total number of
vertical levels
integer :: nlev !Number of vertical levels + sfc
level
integer :: flag !Indicates whether there are
surface variables (=1 if yes)
real :: timeshift !(time of time group)-(time of
this stations obs)
integer :: vvarcounter !Counter over vertical
(multi-level) variables
integer :: vvarcountermax !Total number of vertical
(multi-level) variables
integer :: svarcounter !Counter over surface
(single-level) variables
integer :: svarcountermax !Total number of surface
(single-level) variables
integer :: outunitctl = 13 !Output file unit for
ctl file
character(len=1), parameter :: alpha(1:26) = &
(/'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s',&
&'t','u','v','w','x','y','z'/) ! Alphabet
!Open the output file
open(outunit,file='grads.dat',status='replace',action='write',&
form='unformatted',access='sequential')
!Total number of stations
stationcountermax=3
!Total number of times
timecountermax=5
!Total number of vertical levels in multiple-level
variables
levelcountermax=8
!Total number of surface variables
svarcountermax=2
if(svarcountermax.gt.0) then
!Are there surface variables: =1 if yes, =2 if no
flag=1
else
flag=0
endif
!Number of vertical levels + sfc level
nlev=levelcountermax+flag
!Number of multiple level variables
vvarcountermax=3
if((vvarcountermax.gt.0).and.(levelcountermax.eq.0))
then
print *,'You specified a non-zero number of
vertical variables ',&
'but zero vertical levels.'
print *,'Stopping program due to this error.'
STOP
endif
allocate(lat(stationcountermax))
allocate(lon(stationcountermax))
allocate(stationid(stationcountermax))
!Assign lat/lons and stationids
do stationcounter=1,stationcountermax
lat(stationcounter)=10.0+stationcounter*1.0
lon(stationcounter)=15.0+stationcounter*1.0
write(stationid(stationcounter),'I8'),stationcounter
print *,'stationid
=*',stationid(stationcounter),'*'
enddo
!Loop over times
do timecounter=1,timecountermax
!Loop over stations
do stationcounter=1,stationcountermax
!Station header
write(outunit),stationid(stationcounter),lat(stationcounter),lon(stationcounter),timeshift,nlev,flag
!If there are sfc variables, then write them out
if(flag.eq.1) then
!Sfc data and multi-level data output in same
statement
write(outunit),((2.0**timecounter)*(3.0**stationcounter)*(5.0**svarcounter),svarcounter=1,svarcountermax),&
(levelcounter*1.1,(1.0*levelcounter*vvarcounter,vvarcounter=1,vvarcountermax),&
levelcounter=1,levelcountermax)
else
!Write height of level, followed by the values
of variables at this height
write(outunit),(levelcounter*1.1,&
(1.0*levelcounter*vvarcounter,vvarcounter=1,vvarcountermax),&
levelcounter=1,levelcountermax)
endif
enddo
!Time group terminator
write(outunit),'TIMETERM',0.0,0.0,0.0,0,0
enddo
close(outunit)
print *,'Number of stations = ',stationcountermax
print *,'Number of times = ',timecountermax
print *,'Number of vertical levels =
',levelcountermax
print *,'Are there surface variables = ',flag
print *,'Number of surface variables =
',svarcountermax
print *,'Number of multiple-level variables =
',vvarcountermax
print *,'nlev = ',nlev
!Open ctl
open(outunitctl,file='grads.ctl',status='replace',action='write',&
form='formatted')
write(outunitctl,'A17'),'DSET ^grads.dat'
write(outunitctl,'A14'),'DTYPE station'
write(outunitctl,'A16'),'STNMAP grads.map'
write(outunitctl,'A13'),'UNDEF -999.0'
write(outunitctl,'A16'),'TITLE Test Data'
write(outunitctl,'A18'),'OPTIONS sequential'
write(outunitctl,'A5,I4,A28'),'TDEF
',timecountermax,' linear 00:00Z01may2002 01mn'
write(outunitctl,'A5,I4'),'VARS
',svarcountermax+vvarcountermax
do svarcounter=1,svarcountermax
write(outunitctl,'A1,A20'),alpha(svarcounter),' 0
99 no description'
enddo !Loop over surface variables
do vvarcounter=1,vvarcountermax
write(outunitctl,'A3,A1,A20'),'aaa',alpha(vvarcounter),'
1 99 no description'
enddo !Loop over multiple-level variables
write(outunitctl,'A7'),'ENDVARS'
close(outunitctl)
END PROGRAM TSMIN
__________________________________
Yahoo! Mail - PC Magazine Editors' Choice 2005
http://mail.yahoo.com
More information about the gradsusr
mailing list