program rds4gscnhdf C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Name - rds4gscnhdf.f Type - Program C Version - 1.0 Date - 12/29/92 Programmer - J. Travers, SAIC C Version - 2.0 Date - 05/26/94 Programmer - Langely DAAC C Version - 2.1 Date - 06/29/94 Programmer - Langley DAAC C Version - 2.2 Date - 07/26/94 Programmer - Langely DAAC C C Purpose - This program reads a monthly S-4G scanner HDF file (either 2.5, 5.0, C 10.0 degree regional, or the zonal/global file) C C 05/26/94 The longitude values which are stored in the S-4G HDF files C are being stored from -180.0 to 180.0 where as the data C is being stored from 0 to 360 degrees. C 06/29/94 By using the ichar command on 8-bit values, it only converts C positive numbers to positive numbers correctly. The negative C numbers were also converted to positive values as well. To C correct this problem, an array was assigned in the readdata C subroutine which stores the correct values for the 8-bit C integers that are suppose to be negative values. This problem C only occured in the FORTRAN read program. C C !!!!!!!!!!!!!!!!!!!!!! NOTE !!!!!!!!!!!!!!!!!!!!!!!!!!!! C The longitude values in the data files have not been corrected. C Only this read program has been modified to correct this C presentation of the data. This is only affected when a user C selects data from a latitude/longitude option. C If the user decides to use an outside piece of software which C depends on these latitude/longitude values, please ignore C these values. C C 07/26/94 The minimum and maximum values are incorrect within the HDF data C files for the longwave and shortwave parameters. So that the C user would not get this information from the data file, the C print lines for this information have been commented out in C this code. For the correct values for all parameters, please C refer to the ERBE S-4G User's Guide. C C C Input - HDF file named s4g_sc2.5_yymm_s, s4g_sc5.0_yymm_s, s4g_sc10._yymm_s C s4g_sc_zg_yymm_s C C Output - Report file names are HDF file names as above plus .rpt C C Key Local Parameters : ddate - data date of file C satcd - satellite code of file C hdfname - name of HDF file C outfile - name of report file C numdims - number of dimensions C dimsize - dimension sizes C dataitems - reference numbers of dataitems wanted C type - tells whether user specified subset as C lat and lon or region number or none C begin1- beginning latitude in degrees or region numb C begin2- beginning longitude in degrees or 0 if region C end1 - ending latitude in degrees or region number C end2 - ending longitude in degrees or 0 if region C out - pointer to report file C contkey - continue program or not C datatyp - type of spatial resolution C resol - resolution 2.5, 5.0, 10.0 degrees C C Subroutines Called : ask - ask for input C getfilenm - create input and output file names C readhdf - reads hdf file id and description C getarea - gets area to subset C getdata - gets parameters to subset C readdata - reads data and attributes from HDF file CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC integer maxrank,listsize parameter (maxrank = 2, listsize = 650) integer datatyp, resol integer dataitems(listsize) real begin1, begin2, end1, end2 character type*2, ddate*5, satcd*2, hdfname*50, outfile*50 character contkey write (6,*) write(6,*)'***************************************************' write(6,*)'* S-4G HDF SCANNER READ PROGRAM *' write(6,*)'* *' write(6,*)'* *' write(6,*)'* Version 2.2 July 26, 1994 *' write(6,*)'* *' write(6,*)'* NOTE TO THE USER!!!!! *' write(6,*)'* FOR YOUR INFORMATION AS HOW THE DATA HAS BEEN *' write(6,*)'* STORED IN THE S-4G HDF FILE, LATITUDE VALUES *' write(6,*)'* AND LONGITUDE VALUES HAVE BEEN ADDED TO THE *' write(6,*)'* INPUT FILE. PLEASE IGNORE THESE VALUES WHEN *' write(6,*)'* USING OTHER HDF PACKAGES BECAUSE THESE NUMBERS *' write(6,*)'* WILL NOT REPRESENT THE DATA CORRECTLY. THIS *' write(6,*)'* READ PROGRAM HAS BEEN MODIFIED TO HANDLE THESE *' write(6,*)'* VALUES. IF YOU HAVE ANY QUESTIONS, PLEASE FEEL *' write(6,*)'* FREE TO CONTACT THE LANGLEY DAAC. *' write(6,*)'* THE WAY THAT THE LATITUDE AND LONGITUDE VALUES *' write(6,*)'* ARE STORED IN THE DATA FILE WILL NOT HAVE ANY *' write(6,*)'* MEANING OR REPRESENT THE WAY THAT THE LATITUDE *' write(6,*)'* AND LONGITUDE VALUES ARE STORED IN THE DATA FILE*' write(6,*)'* WILL NOT HAVE ANY MEANING ORE REPRESENT THE DATA*' write(6,*)'* CORRECTLY IN OUTSIDE SOFTWARE PACKAGES. *' write(6,*)'* *' write(6,*)'* 07/26/94 THE MINIMUM AND MAXIMUM VALUES ARE *' write(6,*)'* INCORRECT WITHIN THE HDF DATA FILES FOR THE *' write(6,*)'* LONGWAVE AND SHORTWAVE PARAMETERS. SO THAT THE *' write(6,*)'* USER WOULD NOT GET THIS INFORMATION FROM THE *' write(6,*)'* DATA FILE, THE PRINT LINES FOR THIS INFORMATION *' write(6,*)'* HAVE BEEN COMMENTED OUT IN THIS CODE. FOR THE *' write(6,*)'* CORRECT VALUES FOR ALL PARAMETERS, PLEASE REFER *' write(6,*)'* TO THE ERBE S-4G USER"S GUIDE. *' write(6,*)'* *' write(6,*)'***************************************************' write(6,*)'* *' write (6,*) 'Please enter Q to exit the program or any other' write (6,*) 'character to continue. ' read(5,'(a)') contkey if (contkey.eq.'Q' .or. contkey.eq.'q') then write(6,*) 'Exiting program ...' stop else call ask(ddate,satcd,datatyp,resol) endif C *** Call subroutine to get file names *** call getfilenm(ddate,satcd,datatyp,resol,hdfname,outfile) C *** call subroutine to read hdf file id, etc *** print *, 'Reading file id and description...' print *, ' ' call readhdf(hdfname) C *** call subroutine to get area to subset *** call getarea(type,begin1,begin2,end1,end2,datatyp) C *** call subroutine to get data items to subset *** call getdata(hdfname,dataitems,datatyp,resol) C *** call subroutine to read data *** print *, ' ' print *, 'Output file name is ',outfile print *, ' ' print *, 'Getting Data, Please wait ...' print *, ' ' call readdata(hdfname,dataitems,type,begin1, 1 begin2,end1,end2,datatyp,resol) end subroutine ask(ddate,satcd,datatyp,resol) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Name - ask Type - subroutine C Version - 1.0 Date - 10/20/92 Programmer - J. Travers, SAIC C 2.0 12/29/92 J. Travers, SAIC C C Purpose - This subroutine ask the user for input C C Output - ddate - the data date C satcd - the satellite code C datatyp - type of spatial average (regional, zonal, global) C resol - type of resolution (2.5, 5.0, 10.0 degrees) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC character ddate*5, satcd*2 integer datatyp, resol print *, ' ' write(6,1) 1 format('Enter data date of the data (yymm). Enter Q to quit ',$) read(*,'(a5)') ddate write(6,*) If (ddate(1:1) .eq. 'Q' .or. ddate(1:1) .eq. 'q') then print *, ' ' print *, 'Exiting Program ... ' stop else print *, ' ' 2 print *, 'Enter satellite code of data to read (1-7).' write(6,3) 3 format('Enter H for help or Q to quit ',$) read(*,'(a2)') satcd write(6,*) if (satcd(1:1).eq.'H' .or. satcd(1:1).eq.'h') then print *, ' ' print *,'===================================================' print *, ' Enter data date (yymm), Example: 8504' print *, ' Enter one of the following satellite codes' print *, ' ' print *, ' 1 - NOAA-9 4 - NOAA-9 and NOAA-10' print *, ' 2 - ERBS 5 - ERBS and NOAA-9' print *, ' 3 - NOAA-10 6 - ERBS and NOAA-10' print *, ' 7 - All 3 satellites' print *,'==================================================' print *, ' ' go to 2 else if (satcd(1:1) .eq. 'Q' .or. satcd(1:1) .eq. 'q') then print *, 'Exiting program ...' stop else print *, ' ' print *, '1) Regional 2) Zonal 3) Global' write(6,4) 4 format('Choose data type from above : ',$) read(*,*) datatyp write(6,*) if (datatyp.ne.3) then print *, ' ' print *, '1) 2.5 degree 2) 5.0 degree 3) 10.0 degree' write(6,5) 5 format('Choose resolution from above : ',$) read(*,*) resol write(6,*) else resol = 1 endif endif endif return end C C subroutine getfilenm(ddate,satcd,datatyp,resol,hdfname,outfile) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Name - getfilenm Type - Subroutine C Version - 1.0 Date - 12/29/92 Programmer - J. Travers, SAIC C C Purpose - This subroutine creates the input and output file names C C Input - ddate - data date of data C satcd - satellite code of data C datatyp - type of spatial average of data C resol - resolution of data C C Output - hdfname - name of the hdf file C outfile - name of report file C C Key Local Parameters - NONE C C Subroutines Called - NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC character ddate*5, satcd*2, hdfname*50, outfile*50 integer datatyp, resol if (datatyp .eq. 1 .and. resol .eq. 1) then hdfname = 's4g_sc2.5_' // ddate(1:4) // '_' // satcd(1:1) else if (datatyp .eq. 1 .and. resol .eq. 2) then hdfname = 's4g_sc5.0_' // ddate(1:4) // '_' // satcd(1:1) else if (datatyp .eq. 1 .and. resol .eq. 3) then hdfname = 's4g_sc10._' // ddate(1:4) // '_' // satcd(1:1) else if (datatyp .ge. 2) then hdfname = 's4g_sc_zg_' // ddate(1:4) // '_' // satcd(1:1) endif outfile = hdfname(1:16) // '.rpt' open(20,file=outfile,form='formatted') print *, ' ' print *, 'Input file is ', hdfname print *, 'Output file is ', outfile print *, ' ' return end C C subroutine readhdf(hdfname) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Name - readhdf Type - Subroutine C Version - 1.0 Date - 10/20/92 Programmer - J. Travers,SAIC C 2.0 12/29/92 J. Travers,SAIC C C Purpose - This subroutine reads the file id and file description C C Input - hdfname - the name of the HDF input file C Output - NONE C C C Key Local Parameters : iret - return code from HDF routines C dfile - pointer to HDF file C id - file id C filedesc - file description C C Functions Called : HDF functions to get file id and file description CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C integer Hopen, DFACC_read,dagfid,dagfds,Hclose integer iret, dfile, i, maxlen, first parameter (DFACC_read=1, first=1,maxlen=625) character hdfname*50,id*96,filedesc*625 do 5, i=1,625 filedesc(i:i) = ' ' 5 continue C *** open HDF file *** dfile = Hopen(hdfname,DFACC_read,0) if (dfile.eq.0) then write(20,*) 'Error on opening HDF file' write(20,*) stop endif C *** get file id *** iret = dagfid(dfile,id,maxlen,first) if (iret.eq.-1) then write(20,*) 'Error on getting file id' write(20,*) stop else write(20,*) 'File ID: ',id write(20,*) endif C *** get file description *** iret = dagfds(dfile,filedesc,maxlen,first) if (iret.eq.-1) then write(20,*) 'Error in getting file description' write(20,*) stop else write(20,*) 'File description is as follows:' write(20,*) c j = 1 c do while (j .le. 625) c start = j c write(20, 10) filedesc(j:j) c if (mod(j, 80) .eq. 0 .and. filedesc(j:j) .ne. ' ') then c do while (filedesc(j+1:j+1) .ne. ' ' .and. j .le. 625) c write(20, 10) filedesc(j:j) c j = j + 1 c end do c write(20, *) c endif c filedesc(j:j) = char(10) c j = j + 1 c end do c write(20, *) filedesc(1:625) c endif i = 1 do while (i .le. 625) write(20,10) filedesc(i:i) 10 format(A, $) if (mod(i,80).eq.0) then do while (filedesc(i+1:i+1).ne.' ' .and. i .le. 625) i = i + 1 write(20,10) filedesc(i:i) end do write(20,*) if (filedesc(i+1:i+1) .eq. ' ') then do while(filedesc(i+1:i+1) .eq. ' ' .and. i .le. 625) i = i + 1 end do endif endif i = i + 1 end do endif c *** close file *** iret = Hclose(dfile) return end subroutine getdata(hdfname,dataitems,datatyp,resol) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Name - getdata Type - Subroutine C Version - 1.0 Date - 10/23/92 Programmer - J. Travers,SAIC C 2.0 12/29/92 J. Travers,SAIC C C Purpose - This subroutine allows the user to choose from the annotation C labels, which data parameters to dump C C Input - hdfname - the number of HDF file C datatyp - type of spatial resolution C resol - type of resolution 2.5, 5.0, or 10.0 degrees C Output - dataitems - list of reference numbers of data to dump C C Key Local Parameters - nlabels - number of reference labels C i,j,k,count - counters C itemnum - number for user to choose C startpos - starting position in HDF file C reflist - list of reference numbers C labellist - list of annotation labels C C Functions called - annotation label routines from HDF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C integer nlabels, i, j, k, count, itemnum, dallist, DFTAG_NDG integer maxlng,startpos,listsize,datatyp,resol parameter (maxlng=21, DFTAG_NDG=720,listsize=650) integer dataitems(listsize), reflist(listsize) character hdfname*50, labellist*21 dimension labellist(listsize) C *** initialize dataitems array *** do 40, i=1,listsize dataitems(i) = 65535 40 continue C *** get annotation labels for appropriate data *** if (datatyp.eq.1) then startpos = 1 else if (datatyp.eq.2 .and. resol.eq.1) then startpos = 1 else if (datatyp.eq.2 .and. resol.eq.2) then startpos = 651 else if (datatyp.eq.2 .and. resol.eq.3) then startpos = 1301 else if (datatyp.eq.3) then startpos = 1951 endif nlabels = dallist(hdfname,DFTAG_NDG,reflist,labellist, 1 listsize,maxlng,startpos) if (nlabels.eq.-1) then write(20,*) 'Error in getting annotation labels' write(20,*) stop endif print *, ' ' j = 1 k = 1 i = 1 do while (i .ge. 1 .and. i .le. 650) if (j.eq.1) then itemnum = 1 print *, ' Item Number Label' print *, ' ' endif write(*,41) i,labellist(i) 41 format(15x,i3,14x,a21) C *** Only print one daily label and then skip to next parameter *** if (i.eq.9 .or. i.eq.100 .or. i.eq.131 .or. i.eq.162 .or. 1 i.eq.193 .or. i.eq.224 .or. i.eq.255 .or. i.eq.430 .or. 2 i.eq.461 .or. i.eq.492 .or. i.eq.523) then i = i + 30 endif C *** Only print one hourly label and then skip to next parameter *** if (i.eq.40 .or. i.eq.64 .or. i.eq.286 .or. i.eq.310 .or. 1 i.eq.334 .or. i.eq.376 .or. i.eq.382 .or. i.eq.406 .or. 2 i.eq.554 .or. i.eq.578 .or. i.eq.602 .or. 3 i.eq.626 .or. i.eq.358) then i = i + 23 endif C *** Let user choose parameters *** j = j + 1 if (j.gt.15 .or. i.eq.nlabels) then print *, ' ' print *, 'Please choose the numbers of the data items you wish' print *, 'to dump,one per line. Note that the item numbers' print *, 'skip on daily and hourly parameters to show the ' print *, 'number of days or hours you will receive. Enter ' print *, '-1 to continue to next list or 0 to quit entering' print *, 'parameters.' print *, ' ' do while (itemnum.ne.-1 .and. itemnum.ne.0) read(*,*) itemnum if (itemnum.eq.-1) go to 45 if (itemnum.eq.0) then dataitems(k) = 0 return endif C *** get daily parameter reference numbers *** if (itemnum.eq.9 .or. itemnum.eq.100 .or. itemnum.eq.131 1 .or. itemnum.eq.162 .or. itemnum.eq.193 .or. itemnum.eq.224 2 .or. itemnum.eq.255 .or. itemnum.eq.430 .or. itemnum.eq.461 3 .or. itemnum.eq.492 .or. itemnum.eq.523) then do 43, count = itemnum, itemnum+30 dataitems(k) = reflist(count) k = k + 1 43 continue C *** get hourly parameter reference numbers *** else if (itemnum.eq.40 .or. itemnum.eq.64 .or. itemnum.eq.286 1 .or. itemnum.eq.310 .or. itemnum.eq.334 .or. itemnum.eq.376 2 .or. itemnum.eq.382 .or. itemnum.eq.406 .or. itemnum.eq.554 3 .or. itemnum.eq.578 .or. itemnum.eq.602 .or. 4 itemnum.eq.626) then do 44, count = itemnum, itemnum+23 dataitems(k) = reflist(count) k = k + 1 44 continue else dataitems(k) = reflist(itemnum) k = k+1 endif 45 continue end do j = 1 endif i = i + 1 end do dataitems(k) = 0 return end subroutine getarea(type,begin1,begin2,end1,end2,datatyp) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Name - getarea type - subroutine C Version - 1.0 Date -10/26/92 Programmer - J. Travers,SAIC C 2.0 12/29/92 J. Travers,SAIC C C Purpose - This functions asks the user to enter the latitude and C longitude ranges, region numbers, or band numbers C to subset data C C Input Parameters - datatyp - type of spatial average C C Output Parameters -type - enter lat & lon, regions, or none C begin1-beginning latitude or region number C begin2-beginning longitude or 0 if region C end1 - ending latitude or region number C end2 - ending longitude or 0 if region C C Key Local Parameters -NONE C C Functions Called - NONE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC character type*1 real begin1,begin2,end1,end2 integer datatyp if (datatyp.eq.1) then print *, ' ' print *, 'Enter whether you would like to further subset ' print *, 'based on latitude and longitude (L), region number ' print *, '(R), or no subsetting (N)' read (*,'(a1)') type print *, ' ' if (type.eq.'L' .or. type.eq.'l') then write(*,55) 55 format('Enter beginning and ending latitude (-90 to 90) ',$) read(*,*) begin1,end1 write(*,56) C 05-26-94 longitude values should go from 0 to 360 degrees C 56 format('Enter beginning and ending longitude(-180 to 180) ',$) 56 format('Enter beginning and ending longitude(0 to 360) ',$) read(*,*) begin2,end2 else if (type.eq.'R'.or.type.eq.'r') then write(*,57) 57 format('Enter beginning and ending region number ') write(*,58) 58 format('(1-10368: 2.5 degree) (1-2592: 5.0 degree)', 1 ' (1-648: degree) ',$) read (*,*) begin1,end1 begin2 = 0.0 end2 = 0.0 else if (type.eq.'N' .or. type.eq.'n') then begin1 = 0.0 begin2 = 0.0 end1 = 0.0 end2 = 0.0 else write(*,*) 'Wrong subsetting entry ... Exiting program' stop endif else if (datatyp.eq.2) then write(*,*) 'Enter whether you would like to further subset' write(*,*) 'by latitude band number (B) or no subsetting(N)' read(*,'(a1)') type print *, ' ' if (type.eq.'B' .or. type.eq.'b') then write(*,*) 'Enter beginning and ending band number' write(*,59) 59 format('(1-72: 2.5 degree) (1-36: 5.0 degree)', 1 ' (1-18: 10 degree)',$) read(*,*) begin1, end1 begin2 = 0.0 end2 = 0.0 else if (type.eq.'N' .or. type.eq.'n') then begin1 = 0.0 begin2 = 0.0 end1 = 0.0 end2 = 0.0 else print *, ' ' print *, 'Wrong subsetting enty ... Exiting program' stop endif else type = 'G' end1 = 0.0 end2 = 0.0 begin1 = 0.0 begin2 = 0.0 endif return end subroutine readdata(hdfname,dataitems,type, o begin1,begin2,end1,end2,datatyp,resol) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Name - readdata type - C Version - 1.0 Date-10/26/92 Programmer- J. Travers,SAIC C 2.0 01/04/93 J. Travers,SAIC C 3.0 06/29/94 Langley DAAC C This subroutine has been modified to handle the correct C output of 8-bit values on all platforms. This will also C run on a Cray. C C Purpose - This function gets dimension of SDS, dimension scales and C attibutes, determines the data slice to retrieve and using C the reference numbers gets the appropriate data, data C attributes, max, min, and scale factor. It also applies the C scale factor to the data before printing it out. C C Input Parameters - hdfname - name of hdf file to read C dataitems- stores reference numbers of data chosen C type-type of subsetting (lat,lon, region # or none) C begin1- beginning latitude or region number C begin2- beginning longitude or 0 if region C end1- ending latitude or region number C end2- ending longitude or 0 if region C datatyp - type of spatial average C resol - type of resolution C C Output Parameters - none C C Key Local Parameters - winst - stores starting dimensions of data slice C windims- stores length of slice C dims - stores dimension of data array C data - stores data gotten from HDF file C scale - stores scale factor C datalabel - stores label for data C dataunit- stores units for data C datafmt - stores format for data C sclabel - stores scale factor as annotation text C max - stores maximum value C min - stores minimum value C reg1, reg2 - actual regions obtained in slice C dimsize - size of each dimension C numdims - number of dimensions in SDS C maxrank - max number of dimensions in SDS C dimscal32,16,8 - arrays stores dimension scales C dimlab - dimension label C dimunt - dimension unit C dimfmt - dimension format C numtype - stores datatype of SDS C index - index into dataarray C resoltn - stores resolution 2.5, 5.0, 10.0 C factor - stores factor for dimension scales C iloop1,iloop2 - counters for datarray loop print C latindex - latitude index for output C outindex -index for regions or band output C C C Subroutines Called - HDF routines to read the SDS dimension, dimension C attribute and scales, data slice, max,min C data attributes, etc. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C character hdfname*50, type*1, datalabel*46, datafmt*6 character dataunt*16, datacoor*2, sclabel*5 character dimlab*40, dimunt*30, dimfmt*5 integer maxrank, listsize, dsrref, dsgrang, dsgdast integer dsgnt, dsgslc, dsgdims, daglab,DFTAG_SD,maxlen integer dsgdist, dsgdisc parameter (maxrank=2, listsize=650, DFTAG_SD=702, maxlen=5) integer dimsize(maxrank), dataitems(listsize), winst(maxrank) integer winend(maxrank), windims(maxrank), dims(maxrank) integer scale,i,j,k,iret,numdims,reg1,reg2 integer datatyp, resol, index, numtype, iloop1, iloop2, outindex integer*4 data32(10368), max32, min32, dimscal32(144) integer*2 data16(10368), max16, min16, dimscal16(144) character chdata8(10368), chmax8, chmin8, chdims8(144) integer data8(10368), count1, dimscal8(144) real begin1,begin2, end1, end2,datavalue,factor,resoltn,latindex dimension dimlab(2) integer addres8(256), kk k = 0 do while (k .le. 127) addres8(k+1) = k k = k+1 enddo k = 129 kk = -128 do while (k .le. 255) addres8(k) = kk + 1 k = k + 1 kk = kk + 1 enddo C *** according to reference numbers go to next data set *** k = 1 do while (dataitems(k).ne.0) iret = dsrref(hdfname,dataitems(k)) if (iret.eq.-1) then write(20,*) 'Error on going to specified reference', 1 ' number - ',dataitems(k) write(20,*) stop endif C *** get dimensions of the SDS *** if (k.eq.1) then iret = dsgdims(hdfname,numdims,dimsize,maxrank) if (iret.ne.0) then write(20,*) 'Error on getting dimensions' write(20,*) stop else write(20,*) write(20,*) write(20,*) 'Dimensions are',dimsize(1), ' by',dimsize(2) write(20,*) endif C ** get number type *** iret = dsgnt(numtype) C *** get dimension attributes *** do 10, i=1,numdims iret = dsgdist(i,dimlab(i),dimunt,dimfmt) if (iret.ne.0) then write(20,*) 'Error in getting dimension attributes' write(20,*) stop else write(20,*) 'Description : ',dimlab(i) write(20,*) 'Units : ',dimunt write(20,*) 'Format : ',dimfmt write(20,*) endif 10 continue C *** get dimension scales for all but global *** if (datatyp.ne.3) then do 46, i=1,numdims if (numtype.eq.24) then iret = dsgdisc(i,dimsize(i),dimscal32) else if (numtype.eq.22) then iret = dsgdisc(i,dimsize(i),dimscal16) else if (numtype.eq.20 .and. datatyp.eq.2 .and. 1 resol.eq.3) then iret = dsgdisc(i,dimsize(i),dimscal16) else if (numtype.eq.20) then iret = dsgdisc(i,dimsize(i),chdims8) do 11 count1 = 1,144 dimscal8(count1) = ichar(chdims8(count1)) if (dimscal8(count1) .ge. 129) then dimscal8(count1) = addres8(dimscal8(count1)) endif write(*,*) 'dimscal8 = ',dimscal8(count1) 11 continue endif if (iret.ne.0) then write(20,*) 'Error on getting dimension scales' stop else if (datatyp.eq.1 .and. resol.eq.1) then factor = 2.5 else if (datatyp.eq.1 .and. resol.eq.2) then factor = 5.0 else if (datatyp.eq.1 .and. resol.eq.3) then factor = 10.0 else factor = 1.0 endif write(20,*) 'Dimension number ', i, ' scales: ' write(20,*) do 44, j=1,dimsize(i) if (numtype.eq.24) then if ( i .eq. 1) then write(20,42) ((dimscal32(j)*factor)+180.0) else write(20,42) (dimscal32(j)*factor) endif else if (numtype.eq.22) then if ( i .eq. 1) then write(20,42) ((dimscal16(j)*factor)+180.0) else write(20,42) (dimscal16(j)*factor) endif else if (numtype.eq.20 .and. datatyp.eq.2 .and. 1 resol.eq.3) then if (j .lt. 10 .and. i.eq.2) then if ( i .eq. 1) then write(20,42) ((dimscal16(j)*factor)+180.0) else write(20,42) (dimscal16(j)*factor) endif else if ( i .eq. 1) then write(20,42) ((dimscal16(j)*factor + j)+180.0) else write(20,42) (dimscal16(j)*factor + j) endif endif else if (numtype.eq.20) then if ( i .eq. 1) then write(20,42) int (dimscal8(j)) * factor + 180.0 else write(20,42) int (dimscal8(j)) * factor endif endif 42 format(f8.1,$) if (mod(j,10).eq.0) then write(20,*) endif 44 continue endif write(20,*) write(20,*) 46 continue endif else C *** get dimensions *** iret = dsgdims(hdfname,numdims,dimsize,maxrank) if (iret.ne.0) then write(20,*) 'Error on getting dimensions' write(20,*) stop endif C *** get number type of SDS *** iret = dsgnt(numtype) endif C *** subset by region numbers *** if (type.eq.'R' .or. type.eq.'r') then winst(2) = int (begin1/dimsize(1)) + 1 winst(1) = 1 if (mod(int (end1), dimsize(1)) .eq. 0) then winend(2) = int (end1/dimsize(1)) else winend(2) = int ((end1/dimsize(1)) +1) endif winend(1) = dimsize(1) windims(1) = winend(1) - winst(1) + 1 windims(2) = winend(2) - winst(2) + 1 if (windims(1).gt.dimsize(1)) windims(1) = dimsize(1) if (windims(2).gt.dimsize(2)) windims(2) = dimsize(2) dims(1) = windims(1) dims(2) = windims(2) reg1 = (winst(2) - 1)*dimsize(1) + 1 reg2 = (winend(2)*dimsize(1)) C *** subset by latitude and longitude *** elseif ((type.eq.'L') .or. (type.eq.'l')) then if (resol.eq.1) then resoltn = 2.5 else if (resol.eq.2) then resoltn = 5.0 else if (resol.eq.3) then resoltn = 10.0 endif winst(2) = int (((90.0 - end1)/resoltn) +1) C 05-26-94 to correct longitude to be presented in 0 to 360 C degree window. C winst(1) = int (((-180.0 - begin2)/resoltn)*(-1) + 1) winst(1) = int ((begin2/resoltn) + 1) winend(2) = int (((90.0 - begin1)/resoltn) + 1) C 05-26-94 to correct longitude to be presented in 0 to 360 C degree window. winend(1) = int ((end2/resoltn) + 1) if (winst(2).gt.dimsize(2)) winst(2) = dimsize(2) if (winst(1).gt.dimsize(1)) winst(1) = dimsize(1) if (winend(2).gt.dimsize(2)) winend(2) = dimsize(2) if (winend(1).gt.dimsize(1)) winend(1) = dimsize(1) windims(1) = winend(1) - winst(1) + 1 windims(2) = winend(2) - winst(2) + 1 dims(1) = windims(1) dims(2) = windims(2) C *** no subsetting required or global*** else if (type.eq.'N'.or.type.eq.'n'.or.type.eq.'G' 1 .or.type.eq.'g') then winst(1) = 1 winst(2) = 1 windims(1) = dimsize(1) windims(2) = dimsize(2) dims(1) = windims(1) dims(2) = windims(2) C *** subset by band *** else if (type.eq.'B'.or.type.eq.'b') then winst(2) = int (begin1) winst(1) = 1 winend(1) = 1 winend(2) = int (end1) windims(1) = 1 windims(2) = int (end1 - begin1 + 1) dims(1) = windims(1) dims(2) = windims(2) endif C *** get data attributes *** iret = dsgdast(datalabel,dataunt,datafmt,datacoor) if (iret.ne.0) then write(20,*) 'Error on getting data attributes' write(20,*) stop else write(20,*) 'Data attributes are as follows: ' write(20,*) ' ' write(20,*) 'Label : ',datalabel write(20,*) 'Units : ',dataunt write(20,*) 'Format: ',datafmt write(20,*) ' ' endif C *** get data max and min values for 32 bit data *** if (numtype.eq.24) then iret = dsgrang(max32,min32) if (iret.eq.-1) then write(20,*) 'Error on getting data max and min' write(20,*) stop C 07/26/94 The minimum and maximum values are incorrect within the HDF data C files for the longwave and shortwave parameters. So that the C user would not get this information from the data file, the C print lines for this information have been commented out in C this code. For the correct values for all parameters, please C refer to the ERBE S-4G User's Guide. C else C write(20,*) 'Data Max: ',max32 C write(20,*) 'Data Min: ',min32 C write(20,*) endif C *** get data slice *** iret = dsgslc(hdfname,winst,windims,data32,dims) if (iret.eq.-1) then write(20,*) 'Error on getting data slice' write(20,*) stop else iret = daglab(hdfname,DFTAG_SD,dataitems(k),sclabel,maxlen) if (iret.eq.-1) then write(20,*) 'Error on getting scale factor annotation' write(20,*) stop else read(sclabel,'(i3)') scale write(20,*) 'Scale factor of ', scale, ' has been applied', 1 ' to the data' write(20,*) endif C *** print out data after applying scale factor *** if (type.eq.'R' .or. type.eq.'r') then write(20,48) 'Slice encompassing regions desired starts', 1 ' at region ',reg1,' and ends at region ',reg2 48 format(a,a,i5,a,i5) write(20,*) 'Values are printed by ascending region ', 1 'number from left to right' write(20,*) write(20,*) 'Region Data ' write(20,*) outindex = reg1 else if (type.eq.'L' .or. type.eq.'l') then write(20,50) 'Slice encompasses latitude ',end1,' to', 1 begin1,' and longitude ',begin2, ' to ',end2 50 format(a,f6.1,a,f6.1,a,f6.1,a,f6.1) write(20,*) 'Values are printed in sections representing', 1 ' all requested longitudes for a given','latitude' write(20,*) write(20,*) 'Please note that lat and lon numbers are', 1 ' for the top and left lat and lon of the','cell' write(20,*) write(20,51) 'Therefore the last latitude will be',factor, 1 ' degrees less than what you specified' write(20,*) 51 format(a,f4.1,a) else if (type.eq.'B' .or. type.eq.'b') then write(20,52) 'Slice encompasses zonal values for latitude', 1 ' band ', int (begin1),' to ', int (end1) 52 format(a,a,i2,a,i2) write(20,*) 'Values are printed by ascending band numbers' write(20,*) write(20,*) 'Band # Data' write(20,*) outindex = int(begin1) else if (type.eq.'G' .or. type.eq.'g') then write(20,*) 'Slice is all 3 global values for 2.5, 5.0', 1 ' and 10.0 degrees, respectively' write(20,*) else if ((type.eq.'N'.or.type.eq.'n').and.datatyp.eq.1) then write(20,*) 'All data values are printed out for ', 1 'resolution chosen.' write(20,*) 'Each section represents all longitude', 1 ' values for a given latitude' write(20,*) write(20,*) 'Please note that lat and lon numbers are ', 1 'for the top and left lat and lon of the','cell.' write(20,*) write(20,51) 'Therefore the last latitude will be ', factor, 1 ' degrees less than what you specified' write(20,*) else if ((type.eq.'N'.or.type.eq.'n').and.datatyp.eq.2) then write(20,*) 'All data values are printed out for resolution', 1 ' chosen.' write(20,*) 'Each section represents all latitude bands' write(20,*) write(20,*) 'Band # Data' write(20,*) outindex = 1 endif write(20,*) if (dims(1).eq.1) then iloop1 = dims(1) iloop2 = dims(2) else iloop1 = dims(2) iloop2 = dims(1) endif do 60, i=1, iloop1 write(20,*) if (type.eq.'R' .or. type.eq.'r' .or. type.eq.'B' .or. 1 type.eq.'b' .or. ((type.eq.'N' .or. type.eq.'n') 2 .and. datatyp.eq.2)) then write(20,53) outindex 53 format(i5,$) else if (type.ne.'G' .and. type.ne.'g') then if (i.eq.1) then if (type.eq.'N'.or. type.eq.'n') then latindex = 90.0 C 05-26-94 To correct longitude values to be presented in C 0 to 360 degree format. C begin2 = -180.0 C end2 = 180.0 begin2 = 0.0 end2 = 360.0 else latindex = end1 endif endif write(20,*) 'Data Values are for longitude',begin2, 1 ' to', end2, ' and latitude',latindex write(20,*) latindex = latindex - factor endif do 55, j=1, iloop2 index = (i-1)*iloop2 + j if (data32(index).ne.32767 .and. data32(index) 1 .ne.-32767 .and. data32(index).ne.-127 2 .and. data32(index).ne.2147483647) then datavalue = (float(data32(index))/float(scale)) write(20,54) datavalue 54 format(f17.4,$) else write(20,54) float (data32(index)) endif if (type.eq.'R'.or.type.eq.'r' .or. type.eq.'B' 1 .or.type.eq.'b' .or. ((type.eq.'N'.or.type.eq.'n') 2 .and. datatyp.eq.2)) then outindex = outindex + 1 endif if ((mod(j,5)).eq.0) then write(20,*) if (type.eq.'R'.or.type.eq.'r' .or. type.eq.'B'.or. 1 type.eq.'b' .or. ((type.eq.'N'.or.type.eq.'n') 2 .and. datatyp.eq.2)) then write(20,53) outindex endif endif 55 continue write(20,*) write(20,*) 60 continue endif C *** get data max and min values for 16 bit data *** else if (numtype.eq.22) then iret = dsgrang(max16,min16) if (iret.eq.-1) then write(20,*) 'Error on getting data max and min' write(20,*) stop C 07/26/94 The minimum and maximum values are incorrect within the HDF data C files for the longwave and shortwave parameters. So that the C user would not get this information from the data file, the C print lines for this information have been commented out in C this code. For the correct values for all parameters, please C refer to the ERBE S-4G User's Guide. C else C write(20,*) 'Data Max: ',max16 C write(20,*) 'Data Min: ',min16 C write(20,*) endif C *** get data slice *** iret = dsgslc(hdfname,winst,windims,data16,dims) if (iret.eq.-1) then write(20,*) 'Error on getting data slice' write(20,*) stop else iret = daglab(hdfname,DFTAG_SD,dataitems(k),sclabel,maxlen) if (iret.eq.-1) then write(20,*) 'Error on getting scale factor annotation' write(20,*) stop else if (sclabel(3:3).ne.'0') then read(sclabel,'(i2)') scale else read(sclabel,'(i4)') scale endif write(20,*) 'Scale factor of ', scale,' has been ', 1 'applied to the data' write(20,*) endif C *** print out data after applying scale factor *** if (type.eq.'R' .or. type.eq.'r') then write(20,48) 'Slice encompassing regions desired starts', 1 ' at region ',reg1,' and ends at region ',reg2 write(20,*) 'Values are printed by ascending region ', 1 'number from left to right' write(20,*) write(20,*) 'Region Data ' write(20,*) outindex = reg1 else if (type.eq.'L' .or. type.eq.'l') then write(20,50) 'Slice encompasses latitude ',end1,' to', 1 begin1,' and longitude ',begin2, ' to ',end2 write(20,*) 'Values are printed in sections representing', 1 ' all requested longitudes for a given', 2 'latitude' write(20,*) write(20,*) 'Please note that lat and lon numbers are ', 1 'for the top and left lat and lon of the','cell' write(20,*) write(20,51) 'Therefore the last latitude will be',factor, 1 ' degrees less than what you specified' write(20,*) else if (type.eq.'B' .or. type.eq.'b') then write(20,52) 'Slice encompasses zonal values for latitude', 1 ' band ', int (begin1),' to ', int (end1) write(20,*) 'Values are printed by ascending band numbers' write(20,*) write(20,*) 'Band # Data' write(20,*) outindex = int(begin1) else if (type.eq.'G' .or. type.eq.'g') then write(20,*) 'Slice is all 3 global values for 2.5, 5.0', 1 ' and 10.0 degrees, respectively' write(20,*) else if ((type.eq.'N'.or.type.eq.'n').and.datatyp.eq.1) then write(20,*) 'All data values are printed out for', 1 ' resolution chosen.' write(20,*) 'Each section represents all longitude', 1 ' values for a given latitude' write(20,*) write(20,*) 'Please note that lat and lon numbers are ', 1 'for the top and left lat and lon of the','cell' write(20,*) write(20,51) 'Therefore the last latitude will be ', factor, 1 ' degrees less than what you specified' write(20,*) else if ((type.eq.'N'.or.type.eq.'n').and.datatyp.eq.2) then write(20,*) 'All data values are printed out for', 1 ' resolution chosen.' write(20,*) 'Each section represents all latitude bands' write(20,*) write(20,*) 'Band # Data' write(20,*) outindex = 1 endif write(20,*) if (dims(1).eq.1) then iloop1 = dims(1) iloop2 = dims(2) else iloop1 = dims(2) iloop2 = dims(1) endif do 70, i=1, iloop1 write(20,*) if (type.eq.'R' .or. type.eq.'r' .or. type.eq.'B' .or. 1 type.eq.'b' .or. ((type.eq.'N' .or. type.eq.'n') 2 .and. datatyp.eq.2)) then write(20,53) outindex else if (type.ne.'G' .and. type.ne.'g') then if (i.eq.1) then if (type.eq.'N'.or. type.eq.'n') then latindex = 90.0 C 05-26-94 To present the longitude values in 0 to 360 degree format. C begin2 = -180.0 C end2 = 180.0 begin2 = 0.0 end2 = 360.0 else latindex = end1 endif endif write(20,*) 'Data Values are for longitude',begin2, 1 ' to', end2, ' and latitude',latindex write(20,*) latindex = latindex - factor endif do 65, j=1, iloop2 index = (i-1)*iloop2 + j if (data16(index).ne.32767 .and. data16(index) 1 .ne.-32767 .and. data16(index).ne.-127) then datavalue = (float(data16(index))/float(scale)) write(20,54) datavalue else write(20,54) float (data16(index)) endif if (type.eq.'R'.or.type.eq.'r' .or. type.eq.'B' 1 .or.type.eq.'b' .or. ((type.eq.'N'.or.type.eq.'n') 2 .and. datatyp.eq.2)) then outindex = outindex + 1 endif if ((mod(j,5)).eq.0) then write(20,*) if (type.eq.'R'.or.type.eq.'r' .or. type.eq.'B'.or. 1 type.eq.'b' .or. ((type.eq.'N'.or.type.eq.'n') 2 .and. datatyp.eq.2)) then write(20,53) outindex endif endif 65 continue write(20,*) write(20,*) 70 continue endif C *** get data max and min values for 8 bit data *** else if (numtype.eq.20) then iret = dsgrang(chmax8,chmin8) if (iret.eq.-1) then write(20,*) 'Error on getting data max and min' write(20,*) stop C 07/26/94 The minimum and maximum values are incorrect within the HDF data C files for the longwave and shortwave parameters. So that the C user would not get this information from the data file, the C print lines for this information have been commented out in C this code. For the correct values for all parameters, please C refer to the ERBE S-4G User's Guide. C else c c *** 06-29-94 *** added this check since ichar only works with c positive numbers. If there are negative numbers, now they can c be included in the data output as negative numbers instead of c being converted to positive values. c C if ((ichar(chmax8)) .ge. 129) then C write(20,*) 'Data Max: ',addres8(ichar(chmax8)) C else C write(20,*) 'Data Max: ',ichar(chmax8) C endif C if ((ichar(chmin8)) .ge. 129) then C write(20,*) 'Data Min: ',addres8(ichar(chmin8)) C else C write(20,*) 'Data Min: ',ichar(chmin8) C endif C write(20,*) endif C *** get data slice *** iret = dsgslc(hdfname,winst,windims,chdata8,dims) do 71 count1=1,10368 data8(count1) = ichar(chdata8(count1)) c *** 06-29-94 *** Made changes so that correct output will be generated for the c 8-bit array. ICHAR only generates positive numbers. This test will allow this c software to run on all platforms. if (data8(count1) .ge. 129) then kk = addres8(data8(count1)) else kk = data8(count1) endif data8(count1) = kk 71 continue if (iret.eq.-1) then write(20,*) 'Error on getting data slice' write(20,*) stop else iret = daglab(hdfname,DFTAG_SD,dataitems(k),sclabel,maxlen) if (iret.eq.-1) then write(20,*) 'Error on getting scale factor annotation' write(20,*) stop else read(sclabel,'(i1)') scale write(20,*) 'Scale factor of ', scale,' has been ', 1 'applied to the data' write(20,*) endif C *** print out data after applying scale factor *** if (type.eq.'R' .or. type.eq.'r') then write(20,48) 'Slice encompassing regions desired starts', 1 ' at region ',reg1,' and ends at region ',reg2 write(20,*) 'Values are printed by ascending region', 1 ' number from left to right' write(20,*) write(20,*) 'Region Data ' write(20,*) outindex = reg1 else if (type.eq.'L' .or. type.eq.'l') then write(20,50) 'Slice encompasses latitude ',end1,' to', 1 begin1,' and longitude ',begin2, ' to ',end2 write(20,*) 'Values are printed in sections representing', 1 ' all requested longitudes for a given', 2 'latitude' write(20,*) write(20,*) 'Please note that lat and lon numbers are ', 1 'for the top and left lat and lon of the','cell' write(20,*) write(20,51) 'Therefore the last latitude will be',factor, 1 ' degrees less than what you specified' write(20,*) else if (type.eq.'B' .or. type.eq.'b') then write(20,52) 'Slice encompasses zonal values for latitude', 1 ' band ', int (begin1),' to ', int (end1) write(20,*) 'Values are printed by ascending band numbers' write(20,*) write(20,*) 'Band # Data' write(20,*) outindex = int(begin1) else if (type.eq.'G' .or. type.eq.'g') then write(20,*) 'Slice is all 3 global values for 2.5, 5.0', 1 ' and 10.0 degrees, respectively' write(20,*) else if ((type.eq.'N'.or.type.eq.'n').and.datatyp.eq.1) then write(20,*) 'All data values are printed out for', 1 ' resolution chosen' write(20,*) 'Each section represents all longitude', 1 ' values for a given latitude' write(20,*) write(20,*) 'Please note that lat and lon numbers are ', 1 'for the top and left lat and lon of the','cell' write(20,*) write(20,51) 'Therefore the last latitude will be ', factor, 1 ' degrees less than what you specified' write(20,*) else if ((type.eq.'N'.or.type.eq.'n').and.datatyp.eq.2) then write(20,*) 'All data values are printed out for', 1 ' resolution chosen' write(20,*) 'Each section represents all latitude bands' write(20,*) write(20,*) 'Band # Data' write(20,*) outindex = 1 endif write(20,*) if (dims(1).eq.1) then iloop1 = dims(1) iloop2 = dims(2) else iloop1 = dims(2) iloop2 = dims(1) endif do 80, i=1, iloop1 write(20,*) if (type.eq.'R' .or. type.eq.'r' .or. type.eq.'B' .or. 1 type.eq.'b' .or. ((type.eq.'N' .or. type.eq.'n') 2 .and. datatyp.eq.2)) then write(20,53) outindex else if (type.ne.'G' .and. type.ne.'g') then if (i.eq.1) then if (type.eq.'N'.or. type.eq.'n') then latindex = 90.0 C 05-26-94 To correct the longitude discrepancy and produce C the 0 to 360 degree format. C begin2 = -180.0 C end2 = 180.0 begin2 = 0.0 end2 = 360.0 else latindex = end1 endif endif write(20,*) 'Data Values are for longitude',begin2, 1 ' to', end2, ' and latitude',latindex write(20,*) latindex = latindex - factor endif do 75, j=1, iloop2 index = (i-1)*iloop2 + j if (data8(index).ne.127 .and. data8(index).ne.-127) then datavalue = float (data8(index) / scale) write(20,54) datavalue else write(20,54) float (data8(index)) endif if (type.eq.'R'.or.type.eq.'r' .or. type.eq.'B' 1 .or.type.eq.'b' .or. ((type.eq.'N'.or.type.eq.'n') 2 .and. datatyp.eq.2)) then outindex = outindex + 1 endif if ((mod(j,5)).eq.0) then write(20,*) if (type.eq.'R'.or.type.eq.'r' .or. type.eq.'B'.or. 1 type.eq.'b' .or. ((type.eq.'N'.or.type.eq.'n') 2 .and. datatyp.eq.2)) then write(20,53) outindex endif endif 75 continue write(20,*) write(20,*) 80 continue endif endif k = k +1 end do return end