PROGRAM hdf2grads
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM HDF2GRADS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Program to read history data dump from ARPS HDF 4 format and
! clean-up the dimensions inside the file and create a control file
! for GrADS.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yunheng Wang
! 11/3/2003
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nzsoil ! Grid dimensions.
INTEGER :: nstyps ! Maximum number of soil types.
CHARACTER(LEN=5), PARAMETER :: dim_name(5) = &
(/ 'x_dim','y_dim','z_dim','s_dim','n_dim' /)
CHARACTER(LEN=3), PARAMETER :: mon_name(12) = &
(/ 'jan', 'feb', 'mar', 'apr', 'may', 'jun', &
'jul', 'aug', 'sep', 'oct', 'nov', 'dec' /)
INTEGER :: year, month, day, hour, minute, second
REAL :: time
! INTEGER, PARAMETER :: nhisfile_max=200
! CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
! INTEGER :: hinfmt,nhisfile
! INTEGER :: lengbf,lenfil
! CHARACTER(LEN=80) :: dirname
! NAMELIST /output/ dirname
CHARACTER(LEN=132) :: infilename, outfilename
!-----------------------------------------------------------------------
!
! Misc. Local variables
!
!-----------------------------------------------------------------------
INTEGER :: nf, sindex, aindex
INTEGER :: sd_id, sds_id, dim_id
INTEGER :: nsds, nattr ! File attributes
INTEGER :: srank, stype, snattr ! SD set attributes
INTEGER :: sdim_size(4)
CHARACTER(LEN=20) :: sname
CHARACTER(LEN=20) :: aname
INTEGER :: atype,clen,ulen
CHARACTER(LEN=80) :: comments,units
INTEGER :: ireturn
INTEGER :: n, nvar, nlevel
LOGICAL :: fexist = .TRUE.
INTEGER :: funit, snamelen
CHARACTER(LEN=80) :: vartag(200), timestr
CHARACTER(LEN=20) :: fmtstr
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'hdf.f90'
INTEGER :: sfstart, sfend
INTEGER :: sffinfo, sfginfo
INTEGER :: sfselect, sfendacc
INTEGER :: sfdimid, sfsdmname
INTEGER :: sffattr, sfgainfo, sfrcatt, sfrnatt
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! WRITE(6,'(/9(/5x,a)/)') &
! '###############################################################', &
! '###############################################################', &
! '### ###', &
! '### Welcome to HDF2GRADS ###', &
! '### This program converts the history dump data ###', &
! '### sets generated by ARPS, between various formats. ###', &
! '### ###', &
! '###############################################################', &
! '###############################################################'
!
!-----------------------------------------------------------------------
!
! Get the names of the input data files.
!
!-----------------------------------------------------------------------
!
! CALL get_input_file_names(hinfmt,grdbasfn,hisfile,nhisfile)
WRITE(6,ADVANCE='NO',FMT='(a)') ' Enter HDF file name: '
READ(5,*) infilename
INQUIRE(FILE=trim(infilename), EXIST = fexist )
IF( .NOT. fexist ) THEN
WRITE(6,*) 'ERROR: File "',TRIM(infilename),'" does not exist.'
STOP
END IF
!-----------------------------------------------------------------------
!
! Set the control parameters for output:
!
!-----------------------------------------------------------------------
!
! READ (5,output)
!-----------------------------------------------------------------------
!
! Get dimensions from grid and base file
!
!-----------------------------------------------------------------------
! lengbf = len_trim(grdbasfn)
! CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
! nx,ny,nz,nzsoil,nstyps, ireturn)
!-----------------------------------------------------------------------
!
! Beginning to process HDF files
!
!-----------------------------------------------------------------------
nf = 1
! DO nf = 1,nhisfile ! loop over each file
sd_id = sfstart(infilename,DFACC_WRITE)
IF (sd_id <= 0) THEN
WRITE (6,*) 'ERROR: cannot open file "',TRIM(infilename),'"'
STOP
ENDIF
!
! get dimension information
!
aindex = sffattr(sd_id, 'nx')
ireturn = sfrnatt(sd_id, aindex, nx)
aindex = sffattr(sd_id, 'ny')
ireturn = sfrnatt(sd_id, aindex, ny)
aindex = sffattr(sd_id, 'nz')
ireturn = sfrnatt(sd_id, aindex, nz)
aindex = sffattr(sd_id, 'nzsoil')
ireturn = sfrnatt(sd_id, aindex, nzsoil)
aindex = sffattr(sd_id, 'nstyp')
ireturn = sfrnatt(sd_id, aindex, nstyps)
WRITE(6,'(/a,3(a,I3),2(a,I3))') &
' Dimensions read in from data are:', &
' nx =',nx,', ny =', ny,', nz =',nz, &
' nzsoil =',nzsoil,', nstyps =',nstyps
aindex = sffattr(sd_id, 'year')
ireturn = sfrnatt(sd_id, aindex, year)
aindex = sffattr(sd_id, 'month')
ireturn = sfrnatt(sd_id, aindex, month)
aindex = sffattr(sd_id, 'day')
ireturn = sfrnatt(sd_id, aindex, day)
aindex = sffattr(sd_id, 'hour')
ireturn = sfrnatt(sd_id, aindex, hour)
aindex = sffattr(sd_id, 'minute')
ireturn = sfrnatt(sd_id, aindex, minute)
aindex = sffattr(sd_id, 'second')
ireturn = sfrnatt(sd_id, aindex, second)
aindex = sffattr(sd_id, 'time')
ireturn = sfrnatt(sd_id, aindex, time)
! WRITE(timestr,'(a,I2.2,a1,I2.2,a1,a3,I4.4,F2.0,a2)') &
! '1 LINEAR ',hour,':',minute,'Z', &
! mon_name(month),year,time,'HR'
WRITE(6,'(a,I4.4,a,I2.2,a,I2.2,a,I2.2,a,I2.2,a,I2.2,a,I5.5,a)') &
' Data valid at ', &
year,'-',month,'-',day,'_', &
hour,':',minute,':',second,'+',INT(time),' seconds.'
ireturn = sffinfo(sd_id,nsds,nattr)
nvar = 0
DO sindex = 1,nsds ! loop over each dataset in the file
sds_id = sfselect(sd_id,sindex-1)
ireturn = sfginfo (sds_id,sname,srank,sdim_size,stype,snattr)
!WRITE(6,*) ' ', sname, ' Changed'
!
! 1st dimension
!
dim_id = sfdimid(sds_id,0)
IF(sdim_size(1) == nx) ireturn = sfsdmname(dim_id,dim_name(1))
!
! 2nd dimension
!
dim_id = sfdimid(sds_id,1)
IF(sdim_size(2) == ny) ireturn = sfsdmname(dim_id,dim_name(2))
!
! 3rd dimension
!
dim_id = sfdimid(sds_id,2)
IF(sdim_size(3) == nz) THEN
ireturn = sfsdmname(dim_id,dim_name(3))
ELSE IF(sdim_size(3) == nzsoil) THEN
ireturn = sfsdmname(dim_id,dim_name(4))
ELSE IF(sdim_size(3) == nstyps) THEN
ireturn = sfsdmname(dim_id,dim_name(5))
END IF
!
! 4th dimension
!
dim_id = sfdimid(sds_id,3)
IF(sdim_size(4) == nstyps) THEN
ireturn = sfsdmname(dim_id,dim_name(5))
END IF
!
! Prepare variable tag for GrADS control file
!
aindex = sffattr(sds_id,'comment')
ireturn = sfgainfo(sds_id,aindex,aname,atype,clen)
ireturn = sfrcatt(sds_id,aindex,comments)
aindex = sffattr(sds_id,'units')
ireturn = sfgainfo(sds_id,aindex,aname,atype,ulen)
ireturn = sfrcatt(sds_id,aindex,units)
snamelen = LEN_TRIM(sname)
snamelen = 2*snamelen + 2
WRITE(fmtstr,'(a,I2,a)') '(3a,',22-snamelen,'x,I4,5a)'
IF(srank == 2) THEN
nlevel = 0
IF(sdim_size(1) == nx .AND. sdim_size(2) == ny) THEN
nvar = nvar + 1
WRITE(vartag(nvar),FMT=fmtstr) TRIM(sname),'=>', &
TRIM(sname),nlevel,' 99 ', &
comments(1:clen),' (',units(1:ulen),')'
END IF
ELSE IF(srank == 3) THEN
nlevel = sdim_size(3)
IF(sdim_size(1) == nx .AND. sdim_size(2) == ny .AND. &
sdim_size(3) == nz) THEN
nvar = nvar + 1
WRITE(vartag(nvar),FMT=fmtstr) TRIM(sname),'=>', &
TRIM(sname),nlevel,' 99 ', &
comments(1:clen),' (',units(1:ulen),')'
END IF
ELSE IF(srank == 4) THEN
END IF
ireturn = sfendacc(sds_id)
END DO
ireturn = sfend(sd_id)
!
! Write a GrADS control file
!
WRITE(outfilename,'(a,a)') TRIM(infilename),'.ctl'
WRITE(6,'(/a,a)') ' Output GrADS control file is: ',TRIM(outfilename)
funit = 30 + nf
OPEN(funit, FILE=TRIM(outfilename), FORM='FORMATTED', &
IOSTAT = ireturn, STATUS = 'UNKNOWN')
IF (ireturn /= 0) THEN
WRITE(6,'(1x,3a,I)') 'ERROR: Cannot create file: ', &
TRIM(outfilename),' IOSTAT = ',ireturn
STOP
END IF
WRITE(funit,'(a)') &
'*', &
'* Template of GrADS control file for ARPS HDF format history dumps', &
'*', &
'* Usage:', &
'*', &
'* 1. Use gradshdf instead of grads', &
'* 2. ga-> xdfopen this_control_file', &
'* 3. ga-> set mproj off', &
'* Then works as usual GrADS data files', &
'* 4. Change this control file', &
'* 4.1 Specify dataset: DSET HDF_file_name', &
'* 4.2 Add a variable: VARS current_value + 1', &
'* var_name=>SDS_name Vertical_level 99 comments', &
'* 4.3 Remove a variable: VARS current_value - 1', &
'* remove the variable line', &
'* 5. A perl script in scripts/ directory is recommended.', &
'* 5.1 cd ARPS root directory', &
'* 5.2 source scripts/link_data', &
'* 5.3 cd your working directory', &
'* 5.3 useg -hdf this_control_file', &
'*', &
'* Limitations:', &
'* ', &
'* 1. Does not work with 16bit mapped HDF data', &
'* 2. No map projection, "set mproj off" before start to display any data', &
'* 3. Only works with 2D & 3D atmospheric variables', &
'*', &
'* Note:', &
'* If expects more accurant plots of ARPS data, please use "ARPSCVT" to', &
'* convert HDF files to GrADS binary format.', &
'*', &
'*'
WRITE(funit,'(a,a)') 'DSET ', TRIM(infilename)
WRITE(funit,'(a)' ) 'TITLE GrADS control file for ARPS HDF format dumps'
WRITE(funit,'(a,a,I4,a)') 'XDEF ',dim_name(1), nx,' linear 1 1'
WRITE(funit,'(a,a,I4,a)') 'YDEF ',dim_name(2), ny,' linear 1 1'
WRITE(funit,'(a,a,I4,a)') 'ZDEF ',dim_name(3), nz,' linear 1 1'
! WRITE(funit,'(a,a)') 'TDEF ',TRIM(timestr)
WRITE(funit,'(/a,I4)' ) 'VARS ',nvar
DO n = 1, nvar
WRITE(funit,'(a)') TRIM(vartag(n))
END DO
WRITE(funit,'(a)') 'ENDVARS'
CLOSE(funit)
! END DO
WRITE(6,*)
WRITE(6,'(a)') ' Recommended command sequence to show GrADS plots:',&
' - go to ARPS root directory', &
' - "source scripts/link_data"', &
' - come back to current data directory'
WRITE(6,'(a,a,a/)') ' - "useg -hdf ',TRIM(outfilename),'"'
WRITE(6,*) ' ==== HDF2GRADS terminated normally ===='
END PROGRAM hdf2grads