cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c dxsdsread c c modified by J. Bishop Aug 01 1999 to deal with new name format c and test with HDF 4.3 c and deal with special 8p2 files c modified by J. Bishop Oct 1998 to read all formats and provide c representative ouput. c modified by J. Bishop Jan 02 1996 to read dx sds data sets c c testing program for sds data c modified by J. Bishop Sept 27 1994 c modified by Z. Garraffo Aug 14 1994 c c build with Makefile entry as follows c dxsdsread: dxsdsread.o c f77 dxsdsread.o -L/usr/local/lib -ldf -o dxsdsread c written for HDF 3.3 c c dumps sds header information c and one to two 'image' data sets in ascii format c c ISCCP DX data source data is on an equal area grid with a an c effective 30x30 km pixel size. c For SeaWiFS production, the data have been mapped on a 720 x 360 c (0.5 x 0.5 degree) rectangular grid, with: c c [ 1, 1] centered on 179.75W 89.75N c [ 1,360] centered on " 89.75S c [720, 1] on 179.75E 89.75N c [720,360] on " 89.75S. c c This grid inverts latitude orientation compared to conventional c gcm grid schemes. c The regridding is to facilitate display of sds data sets using c readily available programs such as collage (ncsa) or spyglass. c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc parameter(im=720,jm=360,mrank=2,idebug=1,nspec=3, & rlat=180.,rlng=360.,dlat=rlng/float(im),dlng=rlat/float(jm)) character*80 file,filenam,oldfile character*5 special(nspec),sident character*20 freq character*40 fmt,prefix c 960102 changed to v [character*60 label, units, format, coordsys] character*120 label, units, format, coordsys integer*2 istat integer*4 irank integer dsgdata,dsgdast,dims(2),dsgdims,dsgfill,dsgnt,nt integer*2 isolar(im,jm),max,min,dsgrang real asolar(im,jm) integer*2 fill character fill1,max1,min1,isola1(im,jm) DATA special/'aics8','ctau8','dalb8'/ c c input file name from the command line if it is there c call getarg(1,filenam) if (lnblnk(filenam).lt.1) then write(*,*) ' enter filename ' read(5,101) filenam end if freq = '3 hourly' c ix = index(filenam,'.sds') ix=lnblnk(filenam) ispec=1 do i = 1 , nspec iy=index(filenam,special(i)) if (iy.gt.0) ispec=2 end do write(*,*) 'special =',ispec file=filenam(ix-19:ix) write(*,*) '===========================================' write(*,*) '=' write(*,*) '= File: ',filenam(1:lnblnk(filenam)) 101 format(a80) write(*,*) '=' write(*,*) '= isccp input data version: ',file(1:2) write(*,*) '= data type : ',file(3:6) if (file(7:7).eq.'m') freq = 'monthly' if (file(7:7).eq.'d') freq = 'daily' write(*,*) '= data frequency : ',file(7:7),' => ', & freq write(*,*) '= production version : ',file(8:9) write(*,*) '= year : ',file(11:14) write(*,*) '= month : ',file(15:16) write(*,*) '= HDF data format : ',file(18:20) write(*,*) '=' write(*,*) '= Data should be organized in a ' write(*,*) '= (0.5 x 0.5 degree) rectangular grid:' write(*,*) '=' write(*,*) '= [ 1, 1] centered on 179.75W 89.75N' write(*,80) '= [ 1,',jm,'] centered on " 89.75S' write(*,80) '= [',im,', 1] on 179.75E 89.75N' write(*,85) '= [',im,jm,'] on " 89.75S' write(*,*) '=' write(*,*) '= sds HEADER information follows' igmt = 0 iday = 1 if (file(7:7).eq.'m') iday = 0 if (file(7:7).eq.'8') igmt = 1 c c example: read sds file c istat = DSGDIMS(filenam,irank,dims,mrank) write(*,*) if (istat.ne.0) write(*,*) '= error in dimensions' istat = DSGDAST(label, units, format, coordsys) if (istat.ne.0) write(*,*) '=error getting data specs' write(*,*) '= sds: dims :',dims(1),dims(2) ix=index(label,file(1:lnblnk(file))) if (ix.eq.0) then oldfile=file(1:10)//file(13:20) write(*,*) '=*** label:',label(1:lnblnk(label)) write(*,*) '=*** file name for input:',file(1:lnblnk(file)) ix=index(label,oldfile(1:lnblnk(oldfile))) if(ix.gt.0) then write(*,*) '=*** file name in label :',oldfile(1:18) write(*,*) '=*** WARNING y2k check failed' else write(*,*) '=*** ERROR on ',oldfile(1:lnblnk(oldfile)) stop '=*** ERROR filename does not match label contents' end if end if write(*,90) '= sds: label :',label write(*,90) '= sds: units :',units write(*,90) '= sds: format :',format write(*,90) '= sds:coordsys:',coordsys 80 format(1x,a,i3,a) 85 format(1x,a,i3,',',i3,a) 90 format(1x,a17,a60) ix = index(units,'*') read(units(ix+1:ix+5),91,err=92) imult 91 format(i5) 92 continue if (imult.eq.0) imult=1 write(*,*) '= sds: data are scaled up by ',imult amult = imult istat = dsgnt(nt) write(*,*) '= sds: numtype:',nt if (nt.eq.20) then write(*,*)'= data is signed 8bit integer' else if (nt.eq.22) then write(*,*)'= data is signed 16bit integer' end if istat = dsgfill(fill) write(*,*) '=' write(*,*) '= results data converted to proper scale follow' write(*,*) '=' write(*,*) '= sds: fill :',fill,' => ',float(fill)/amult istat = dsgrang(max,min) write(*,*)'= sds: min,max:',min,max write(*,*)'= sds: => ',float(min)/amult,float(max)/amult write(*,*) '=' mi=999999 do while (istat.eq.0) istat=dsgdata(filenam,2,dims,isolar) if (istat.ne.0) goto 5000 ii = ii+1 c write(*,*)' ii = ',ii,' istat =',istat mx=-999999 mb=999999 do iii = 1,dims(1) do jjj = 1,dims(2) iss=isolar(iii,jjj) mx = max0(mx,iss) mb = min0(mb,iss) if (iss.ge.0) mi= min0(mi,iss) end do end do write(*,*) '= MONTH = ',file(15:16) if (igmt.eq.1) then iiday = int(((ii-1)/ispec)/8)+1 iigmt = mod((ii-1)/ispec,8)*3 write(*,*) '= DAY =',iiday,' GMT = ',iigmt else if (iday.eq.1) then write(*,*) '= DAY =',ii end if if (ii.le.6.and.istat.eq.0) then inum9=dims(1) if (idebug.eq.1) inum9=8 iiint=dims(1)/inum9 write(*,998) ' Latitude y_coord x:',(i8,i8=1,dims(1),iiint) 998 format(a,20i6) 999 format(a,i3,a) if (imult.eq.10000) write(fmt,999)'(a21,',inum9,'(1x,f5.4))' if (imult.eq.1000) write(fmt,999)'(a21,',inum9,'(1x,f5.3))' if (imult.eq.100) write(fmt,999)'(a21,',inum9,'(1x,f5.2))' if (imult.eq.10) write(fmt,999)'(a21,',inum9,'(1x,f5.1))' if (imult.eq.1) write(fmt,999)'(a21,',inum9,'(1x,f5.0))' c write(*,*) fmt, ' = format' do iy= 1,dims(2) alat = -(float(iy)*dlat - dlat/2. - rlat/2.) do i9 = 1,dims(1) asolar(i9,iy) = float(isolar(i9,iy))/float(imult) end do if ((iy-1)/10*10.eq.iy-1) then write(prefix,900) alat,iy 900 format(' ',f6.2,i7,' ') write(*,fmt) prefix,(asolar(ix,iy),ix=1,dims(1),iiint) end if end do end if c end do 5000 write(*,*) '= ' write(*,*) '= data records total found = ',ii write(*,*) '= end of data' write(*,*) '=--------------------------------------------' write(*,*) ' ' stop end integer function lnblnk(msg) character*(*) msg do 5 i=len(msg),1,-1 if (msg(i:i).ne.' ') goto 10 5 continue 10 lnblnk=i return end