Table of Contents 
Previous Chapter
program fortran_sample
c This is a sample fortran program for reading isccp flux dataset
implicit none
integer*4 NCOL,NROW,NMONTH, read_sub
character*256 hdf_file
parameter(NCOL=67,NROW=67, NMONTH=90)
c Meta Data variables
character*256 project
character*256 dataset_name
character*256 source_name
character*256 processing_level
character*256 start_date
character*256 stop_date
character*256 position_type
character*256 grid_type
character*256 grid_name
character*256 projection
real*4 latitude(4)
real*4 longitude(4)
character*256 Temporal_Res
character*256 Spatial_X_Res
character*256 Spatial_Y_Res
c Scientific Data set variables
real*4 latitude_grid(NCOL, NROW)
real*4 longitude_grid(NCOL, NROW)
real*4 dwnvssrf(NCOL, NROW, NMONTH)
real*4 dwnirsrf(NCOL, NROW, NMONTH)
real*4 upvssrf(NCOL, NROW, NMONTH)
real*4 upirsrf(NCOL, NROW, NMONTH)
real*4 dirctop(NCOL, NROW, NMONTH)
real*4 upvstop(NCOL, NROW, NMONTH)
real*4 upirtop(NCOL, NROW, NMONTH)
integer*4 ret
hdf_file ='../hdf/arctic-isccp-fluxes1983-1990.hdf'
c------------------------------------------------------------------------
ret= read_sub(hdf_file,project,dataset_name, source_name,
& processing_level,start_date,
& stop_date,position_type,grid_type,
& grid_name,projection,latitude,longitude,
& Temporal_Res,Spatial_X_Res,Spatial_Y_Res,
& longitude_grid,latitude_grid,dwnvssrf,dwnirsrf,
& upvssrf,upirsrf,dirctop,upvstop,upirtop)
print*,project,len(project)
print*,dataset_name
print*,source_name
print*,processing_level
print*,start_date
print*,stop_date
print*,position_type
print*,grid_type
print*,grid_name
print*,latitude
print*,longitude
print*,Temporal_Res
print*,Spatial_X_Res
print*,Spatial_Y_Res
open(10,file='latgrid.ftst',recl=67*67*4,access='direct',
& form='unformatted')
write(10,rec=1)latitude_grid
close(10)
open(10,file='longrid.ftst',recl=67*67*4,access='direct',
& form='unformatted')
write(10,rec=1)longitude_grid
close(10)
open(10,file='dwnvssrf.ftst',recl=67*67*4*90,access='direct',
& form='unformatted')
write(10,rec=1)dwnvssrf
close(10)
open(10,file='dwnirsrf.ftst',recl=67*67*4*90,access='direct',
& form='unformatted')
write(10,rec=1)dwnirsrf
close(10)
open(10,file='upirsrf.ftst',recl=67*67*4*90,access='direct',
& form='unformatted')
write(10,rec=1)upirsrf
close(10)
open(10,file='upvssrf.ftst',recl=67*67*4*90,access='direct',
& form='unformatted')
write(10,rec=1)upvssrf
close(10)
open(10,file='dirctop.ftst',recl=67*67*4*90,access='direct',
& form='unformatted')
write(10,rec=1)dirctop
close(10)
open(10,file='upirtop.ftst',recl=67*67*4*90,access='direct',
& form='unformatted')
write(10,rec=1)upirtop
close(10)
open(10,file='upvstop.ftst',recl=67*67*4*90,access='direct',
& form='unformatted')
write(10,rec=1)upvstop
close(10)
end
read_sub.f
c---- read_sub----------------------------------------------------
integer function read_sub(hdf_file,project,dataset_name,
& source_name,processing_level,start_date,
& stop_date,position_type,grid_type,
& grid_name,projection,latitude,longitude,
& Temporal_Res,Spatial_X_Res,Spatial_Y_Res,
& longitude_grid,
& latitude_grid,dwnvssrf,dwnirsrf,upvssrf,
& upirsrf,dirctop,upvstop,upirtop)
c---------------------------------------------------------------------------
c This program provides a fortran interface to the TOVS pathP files in HDF
c format.
c
implicit none
integer NROW,NCOL,NMONTH,DFACC_RDONLY
parameter(DFACC_RDONLY=1,NCOL=67,NROW=67, NMONTH=90)
c input:
character *(*) hdf_file
c output:
c Meta Data
character*256 project
character*256 dataset_name
character*256 source_name
character*256 processing_level
character*256 start_date
character*256 stop_date
character*256 position_type
character*256 grid_type
character*256 projection
character*256 grid_name
real*4 latitude(4)
real*4 longitude(4)
character*256 Temporal_Res
character*256 Spatial_X_Res
character*256 Spatial_Y_Res
c scientific data sets
real*4 latitude_grid(NCOL, NROW)
real*4 longitude_grid(NCOL, NROW)
real*4 dwnvssrf(NCOL, NROW, NMONTH)
real*4 dwnirsrf(NCOL, NROW, NMONTH)
real*4 upvssrf(NCOL, NROW, NMONTH)
real*4 upirsrf(NCOL, NROW, NMONTH)
real*4 dirctop(NCOL, NROW, NMONTH)
real*4 upvstop(NCOL, NROW, NMONTH)
real*4 upirtop(NCOL, NROW, NMONTH)
c HDF function declarations
integer*4 sfstart,sffinfo
integer*4 ppgetattr_c,ppgetattr_f,sfend
c internal variables
integer*4 hdf_id,ret,ndatasets,nglobal_attr,error,ret
integer*4 ppgetsds
c Memory storage notes:
c C:
c buffer [NMONTH][NROW][NCOL]
c Fortran,IDL:
c buffer(NCOL,NROW,NMONTH)
c initializations
error=0
c open SD interface to access file with name HDF_FILE
hdf_id =sfstart(HDF_FILE, DFACC_RDONLY)
c You can obtain information about the hdf file through this call
ret = sffinfo(hdf_id,ndatasets, nglobal_attr)
c read Meta data stored as SD attribtutes
ret = ppgetattr_c(hdf_id,"PROJECT", project)
error = error+ret
ret = ppgetattr_c(hdf_id,"DATASET_NAME", dataset_name)
error = error+ret
ret = ppgetattr_c(hdf_id,"SOURCE_NAME", source_name)
error = error+ret
ret = ppgetattr_c(hdf_id,"PROCESSING_LEVEL", processing_level)
error = error+ret
ret = ppgetattr_c(hdf_id,"START_DATE", start_date)
error = error+ret
ret = ppgetattr_c(hdf_id,"STOP_DATE", stop_date)
error = error+ret
ret = ppgetattr_c(hdf_id,"POSITION_TYPE", position_type)
error = error+ret
ret = ppgetattr_c(hdf_id,"GRID_TYPE", grid_type)
error = error+ret
ret = ppgetattr_c(hdf_id,"PROJECTION", projection)
error = error+ret
ret = ppgetattr_c(hdf_id,"Temporal_Res", Temporal_Res)
error = error+ret
ret = ppgetattr_c(hdf_id,"Spatial_X_Res", Spatial_X_Res)
error = error+ret
ret = ppgetattr_c(hdf_id,"Spatial_Y_Res", Spatial_Y_Res)
error = error+ret
ret = ppgetattr_f(hdf_id,"LATITUDE", latitude)
error = error+ret
ret = ppgetattr_f(hdf_id,"LONGITUDE", longitude)
error = error+ret
ret=ppgetsds(hdf_id,"LATITUDE_GRID",latitude_grid,1)
error=error+ret
ret=ppgetsds(hdf_id,"LONGITUDE_GRID",longitude_grid,1)
error=error+ret
ret=ppgetsds(hdf_id,"DWNVSSRF",dwnvssrf,90)
error=error+ret
ret=ppgetsds(hdf_id,"DWNIRSRF",dwnirsrf,90)
error=error+ret
ret=ppgetsds(hdf_id,"DIRCTOP",dirctop,90)
error=error+ret
ret=ppgetsds(hdf_id,"UPVSTOP",upvstop,90)
error=error+ret
ret=ppgetsds(hdf_id,"UPIRTOP",upirtop,90)
error=error+ret
ret=ppgetsds(hdf_id,"UPIRSRF",upirsrf,90)
error=error+ret
ret=ppgetsds(hdf_id,"UPVSSRF",upvssrf,90)
error=error+ret
if (error.lt.0) then
write(0,*) 'Error reading file : ', hdf_file
endif
ret=sfend(hdf_id)
return
end
integer function ppgetsds(hdf_id, sdsname,sds,nz)
c-------------------------------------------------------------------------------
c This function obtains a single scientific data set in file pointed to by
c hdf_id. The name of the scientific data sets is given in sdsname
c
c Input:
c integer*4 hdf_id HDF file ID as returned by sfstart
c character *(*) name of Scientific Data set
c real*4 sds(MAXSIZE)
c This buffer must be large enough to
c hold the entire data set
c integer*4 nz number of elements in z dimension
implicit none
integer*4 MAXSIZE
parameter(MAXSIZE=67*67*90)
integer*4 hdf_id,nz,i,dimsizes(3)
character*(*) sdsname
real*4 sds(MAXSIZE)
c internal variabels
integer*4 index,sds_id,ret,number_type,nattrs,rank
integer*4 START(3),STRIDE(3),EDGE(3)
c hdf function declarations
integer*4 sfn2index,sfrdata,sfselect,sfendacc,sfginfo
c WARNING: STRIDE is set to 1,1,1 due to a bug in HDF3.3r3 the correct
c values should be 0,0,0. Future versions of HDF may have this fixed
data (START(i),i=1,3) /0,0,0/
data (STRIDE(i),i=1,3) /1,1,1/
data (EDGE(i),i=1,3) /67,67,0/
EDGE(3)= nz
c get index of data set in HDF file
index = sfn2index(hdf_id,sdsname)
c select data data set
sds_id = sfselect(hdf_id, index)
c The below call would allow you to obtain information about a
c specific scientific data set
ret = sfginfo(sds_id,sdsname, rank, dimsizes,
& number_type, nattrs)
ppgetsds = sfrdata(sds_id, START, STRIDE, EDGE,sds)
call heprnt(0)
ret = sfendacc(sds_id)
return
end
integer function ppgetattr_c(hdf_id,attribute_name,attribute)
C-----------------------------------------------------------------------
c This function obtains a single global character attributes from the
c PATH P hdf file which is pointed to by hdf_id
c
c input:
c integer*4 hdf_id id returned from sfstart
c character*(*) attribute_name Metadata attribute to return
c character*(*) attribute value of attribute
c returns:
c 0 : success
c -1: error
c-----------------------------------------------------------------------
implicit none
c interface
integer*4 hdf_id
character*(*) attribute_name,attribute
c HDF functions
integer*4 sfrattr,sffattr,sfgainfo
c internal variables
integer*4 index,number_type,length,ret
index=sffattr(hdf_id,attribute_name)
ret = sfgainfo(hdf_id,index,attribute_name, number_type,
& length)
print*,length
ppgetattr_c=sfrattr(hdf_id,index,attribute(1:length))
return
end
integer function ppgetattr_f(hdf_id,attribute_name,attribute)
C-----------------------------------------------------------------------
c This function obtains a 4 element floating point attribute from the
c PATH P hdf file which is pointed to by hdf_id
c
c input:
c integer*4 hdf_id id returned from sfstart
c character*(*) attribute_name Metadata attribute to return
c real*4 attribute(4) value of attribute
c returns:
c 0 : success
c -1: error
c-----------------------------------------------------------------------
implicit none
c interface
integer*4 hdf_id
character*(*) attribute_name
real*4 attribute
c HDF functions
integer*4 sfrattr,sffattr
c internal variables
integer*4 index
index=sffattr(hdf_id,attribute_name)
ppgetattr_f=sfrattr(hdf_id,index,attribute)
return
end
plotgrid_all.pro (IDL routine)
pro plotgrid_all,infile,ngrids
openr,lun,infile,/get_lun
grids=assoc(lun,fltarr(67,67, 90))
month=strarr(12)
months=['January', 'February', 'March', 'April', 'May', 'June', 'July',$
'August', 'September', 'October', 'November', 'December']
pretitle='Fluxes!C'
;levels=findgen(30)*10.
levels=[0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340, $
400,450]
!Y.MARGIN=[3,4]
;set_plot,'ps'
k = 6
year = 1983
FOR i = 0, ngrids-1 DO BEGIN
title = pretitle + months(k) + ' ,' + string(year)
file = strmid(infile, 0, strpos(infile, '.dat'))+'-'$
+strmid(months(k),0,3)+'.ps'
IF (k EQ 11) THEN BEGIN
k = 0
year = year+1
ENDIF ELSE BEGIN
k = k+1
ENDELSE
print, file
; device,file=file,xsize=9,ysize=9,xoffset=5,yoffset=11
FOR n = 0, 7 DO BEGIN
grid3d = grids(i)
grid2d = grid3d(*, *, n)
ntitle = title + 'grid # '+ string(n)
ch, grid2d, ntitle, levels
hak
erase
ENDFOR
erase
; device,/close
end
return
end
pro ch,grid,titlestring,levels
; Functionality : plots a contour plot of values in grid over a polar
; stereographic map. Set up to plot input grid for
; flato-hibler sea ice model
; Input : grid float(36,32) vaues to be contoured
; titlestring string Title string to be printed over the
; plot
; levels(NLEVELS)
; float levels at which contour lines will be drawn
; Maximum number of levels is 30
;
; read input grid
xsize=36 & ysize=32
; levels=[25,50,75,100,125,150,175,200,225,250,275,300,325,350,375,400]
nlevels=n_elements(levels)
annotation=string(format='(F6.0)',levels)
pos=fltarr(4)
latmin_map=60.
rot=0.
!P.CHARSIZE=0.75
; create map
map_set,90.,rot,0.,/stereo,/continents,/grid,limit=[latmin_map,-180.,90.,180.],$
latdel=10,londel=20,mlinestyle=0,mlinethick=5.0,glinestyle=1,$
glinethick=0.5,/noerase,color=128,title=titlestring
;get the map boundaries in normal coordinates.
rot=!map.phioc
lonbottom=rot
lontop=rot+180.
lonleft=rot-90.
lonright=rot+90
cll = fltarr(4, 2)
cll(0, 0:1) = [-90., 60.]
cll(1, 0:1) = [180., 60.]
cll(3, 0:1) = [90., 60.]
cll(2, 0:1) = [0, 60.]
min1=convert_coord(cll(0,0),cll(0,1),/to_norm)
max1=convert_coord(cll(1,0),cll(1,1),/to_norm)
min2=convert_coord(cll(2,0),cll(2,1),/to_norm)
max2=convert_coord(cll(3,0),cll(3,1),/to_norm)
x0=min1(0)
x1=max2(0)
y0=min2(1)
y1=max1(1)
; stop
; set position array for contour overlay
pos = [x0, y0, x1, y1]
print, pos
; perform contouring
; contour,grid,max_val=999.,levels=levels,c_labels=levels,$
; /noerase, $
; position=pos,charsize=1.4, $
; C_linestyle=[0,4],C_charsize=0.8,c_thick=[5,5],/norm,$
; xstyle=4,ystyle=4,c_annotation=annotation
contour,grid,max_val=999,/noerase,position=pos,$
xstyle=5,ystyle=5,/norm,$
C_linestyle=[0,4],C_charsize=0.8,c_thick=[5,5],$
charsize=1.4,levels=levels,c_labels=levels,$
c_annotation=annotation
end
Table of Contents 
Next Chapter