[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: 950711: NetCDF question

>Organization: UC Davis
>Keywords: 199507102313.AA13517 netCDF

Hi Bryan,

> Below is a sample program to write out a netcdf file in VMS. The full data
> set is ECLR(nlon,nlat,ntime). I am trying to write it out as ntime
> matrices of E(nlon,nlat) and then be able to read in the same way. The
> write program has a run time error. I have no trouble writing the full
> data set all at once.  I suspect I don't understand one of the parameters
> in the hyperslab writes and reads.  Thanks for your help

I will need either more information about what the runtime error was or a
stand-alone example program that I can run here to reproduce the error.  I
don't have access to a VMS system, so the latter may be difficult unless you
can reproduce the error on a UNIX system.

>        program wrncdf   !ERBE example
>       PARAMETER (MAXF=40)
>       parameter (nlat=72,nlon=144,ndim=2,ntime=48)

It appears to me that you should be using ndim=3 to define the netCDF
variable 'eclr' with 3 dimensions, nlon, nlat, and ntime, instead of just
the two dimensions nlon and nlat.  Then write one (nlon,nlat) slab for each
time.

>       include 'netcdf.inc'
> C LINK WITH CNVRTC2
> c formats necessary for input indexed files
>       CHARACTER*120 DATAT
>       CHARACTER*45 FILEN(MAXF)
>       CHARACTER*30 DES
>       CHARACTER*10 SHNAME(MAXF)
>       CHARACTER*15 LABX,LABY
>       CHARACTER*8 KEYA
>       CHARACTER*5 UNITS(MAXF),SOURCE(MAXF),AMON(13),LANDSEA(0:2)
>       INTEGER LENGTH(MAXF),IBYTE(MAXF),NVAR(MAXF),NLVL(MAXF)
>       REAL MULT(MAXF)
>       real data(nlon,nlat)
>       DATA IDY/99/,IGM/9/
> c NetCDF definitions
> c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
>       integer  iret
> c* netCDF id
>       integer  ncid
> c* dimension ids
>       integer  latdim, londim ,frtimedim
> c* variable ids
>       integer  Zid, latid, lonid, frtimeid
> c* variable shapes
>       integer dims(ndim)
> c* corners and edge lengths
>       integer corner(ndim), edges(ndim)
> c* data variables ; number of dimensions equal ndim
>       integer*2 idata(nlon,nlat)

This is fine, since you will still only need to keep data for one time in
memory at once.

>       real lat(nlat),lon(nlon),frt
> c* attribute vectors
>       integer*2 intval(2)
> c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
> c GET THE DATA FILES TO BE written
>       OPEN(30,FILE='erbe.LIS',FORM='FORMATTED',STATUS='OLD')
>         READ(30,100)DATAT
>  100    FORMAT(A120)
>         TYPE 100,DATAT
>         N=1
>         NMAPS=1
>  10   READ(30,110,END=15)DES,SHNAME(N),SOURCE(N),UNITS(N),
>      &  LENGTH(N),IBYTE(N),MULT(N),NVAR(N),NLVL(N),FILEN(N)
>  110  FORMAT(A30,A10,2A5,I4,I1,F6.2,I3,I2,A45)
>       TYPE 115, N,DES,SHNAME(N),SOURCE(N),UNITS(N),
>      &  LENGTH(N),IBYTE(N),MULT(N),nvar(n),NLVL(N)
>  115  FORMAT(1x,I2,A30,A10,2A5,I4,I1,F6.2,i3,I2)
>       N=N+1
>       GO TO 10
>  15   CONTINUE
>       CLOSE (30)
> C FIRST READ of ISCCP ERBE DATA FROM FILE 30
>       TYPE *, 'WHICH IS THE ERBE DATA SET YOU WANT?'
>       ACCEPT *, NSELA
>       TYPE *, LENGTH(NSELA)
>       NFa=31
>       OPEN(NFA,FILE=FILEN(NSELA),
>      1  FORM='UNFORMATTED',STATUS='old',ORGANIZATION='INDEXED',
>      1  ACCESS='KEYED',RECL=length(NSELA),KEY=(1:8:CHARACTER))

The above OPEN statement uses VMS-isms that prevent me from being able to
reproduce the error you are seeing here.

> c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
> * enter define mode
>       ncid = nccre ('eclr.dat', NCCLOB, iret)
> * define dimensions
>       latdim = ncddef(ncid, 'lat', nlat, iret)
>       londim = ncddef(ncid, 'lon', nlon, iret)
> c      frtimedim = ncddef(ncid, 'frtime', 3, iret)
>       frtimedim = ncddef(ncid, 'frtime', NCUNLIM, iret)
>       type *, 'finished define dimensions'
> * define variables
>       dims(2) = latdim
>       dims(1) = londim

and also

        dims(3) = frtimedim

>       Zid = ncvdef (ncid, 'eclr', NCshort, ndim, dims, iret)
>       dims(1) = latdim
>       latid = ncvdef (ncid, 'lat', NCFLOAT, 1, dims, iret)
>       dims(1) = londim
>       lonid = ncvdef (ncid, 'lon', NCFLOAT, 1, dims, iret)
>       dims(1) = frtimedim
>       frtimeid = ncvdef (ncid, 'frtime', NCFLOAT, 1, dims, iret)
>       type *, 'finished define variables'
> * assign attributes
>       call ncaptc(ncid, Zid, 'long_name', NCCHAR, 24, 'clear sky long 
>      &wave flux', iret)
>       call ncaptc(ncid, Zid, 'units', NCCHAR, 4, 'W/m2', iret)
>       intval(1) = -500
>       intval(2) = 500
>       call ncapt(ncid, Zid, 'valid_range', NCshort, 2, intval, iret)
>       intval(1) = -999
>       call ncapt(ncid, Zid, '_FillValue', NCshort, 1, intval, iret)
>       call ncaptc(ncid, latid, 'long_name', NCCHAR, 8, 'latitude', iret)
>       call ncaptc(ncid, lonid, 'long_name', NCCHAR, 9, 'longitude', iret
>      1)
>       call ncaptc(ncid, frtimeid, 'long_name', NCCHAR, 10, 'month.year
>      1', iret)
>       type *, 'finished assign attributes'
> * leave define mode
> c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
>       call ncendf(ncid, iret)
> * store lon
>       corner(1) = 1
>       edges(1) = nlon
>       do i=1,nlon
>       ai=i-1
>       lon(i)=21.25 + ai*2.5
>       enddo
>       call ncvpt(ncid, lonid, corner, edges, lon, iret)
>       type *, 'finished store longitude'
> * store lat
>       corner(1) = 1
>       edges(1) = nlat
>       do j=1,nlat
>       aj=j-1
>       lat(j)=-88.75+aj*2.5
>       enddo
>       call ncvpt(ncid, latid, corner, edges, lat, iret)
>       type *, 'finished store latitude'
> c read in data from direct access file
>       itime=0
> c write out main data set
> c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
> c set area of write variables
>       corner(1) = 1
>       corner(2) = 1

and you also need to set corner(3) to specify which time slab to write,
but you will have to do that in the loop, after itime is incremented.


>       edges(1) = nlon
>       edges(2) = nlat

and
        edges(3) = 1

since you only one to write out one time value at a time.

>       DO IYR=85,88
>       DO IMO=1,12
>       type *, 'started store data record'
>       itime=itime+1

        corner(3) = itime

> * store frtime
>       frt=float(imo)+.01*float(iyr)
>       call ncvpt1(ncid, frtimeid, itime, frt,iret)
>       CALL READ25(NSELA,NFA,IBYTE(NSELA),NVAR(NSELA),MULT(NSELA),
>      &  NLVL(NSELA),DATA,IYR,IMO,NLAT,NLON) 
>       do i=1,nlon
>       do j=1,nlat
>       idata(i,j)=data(i,j)
>       enddo
>       enddo
>       call ncvpt(ncid, Zid, corner, edges,idata, iret)
>       type *, 'finished store data record'
>       enddo   !month
>       enddo   !year
>       call ncclos (ncid, iret)
> c $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
>       stop
>       end
>       SUBROUTINE READ25(NSEL,NF,IBYTE,NVAR,MULT,ilev,
>      &  DATA,IYEAR,IMONTH,NLAT,NLON) 
>       PARAMETER(NLO=144,NLA=72)
>       REAL DATA(NLON,NLAT),MULT
>       INTEGER*2 IDUM
>       BYTE B(2)
>       CHARACTER*8 KEYA
>       EQUIVALENCE (IDUM,B(1))
>       DATA IDY/99/,IGM/9/
> C INDEX STRUCTURE
>       STRUCTURE/INPUT2/
>         CHARACTER*8 header
>         integer*2 idat(nlo,nla)
>       END STRUCTURE
>       RECORD/INPUT2/IDAT2
>       STRUCTURE/INPUT1/
>         CHARACTER*8 header
>         CHARACTER*1 idat(nlo,nla)
>       END STRUCTURE
>       RECORD/INPUT1/IDAT1
> C READ DATA
> C SET UP INDEX KEY
>         write(keyA,200) iyear,imonth,idy,igm,ilev
>  200  format(3i2,2i1)
> C        type 210, keyA
>  210  format(2x,a8)
> C READ THE APPROPRIATE DATA
>       IF(IBYTE.EQ.1)THEN
>       IF(ILEV.LT.9)THEN
> C FOR ECMWF UNPACK AND CONVERT
>       read(NF,key=keyA,keyid=0) idat1
>       DO I=1,NLON
>       DO J=1,NLAT
>       IDUM=ICHAR(IDAT1.IDAT(I,J))
>       data(i,j)=b(1)
>       DATA(I,J)=(DATA(I,J) + 128.)*MULT
>       if(data(i,j).lt.0)data(i,j)=-999.
>       ENDDO
>       ENDDO
>       ELSE
> C FOR ISCCP MUST UNPACK AND CONVERT
>       read(NF,key=keyA,keyid=0) idat1
> C CONVERT AND  MULTIPLY BY FACTOR
> c     CALL CONVRTC2(DATA,IDAT1.IDAT,NVAR)
>         do i=1,NLON
>       do j=1,NLAT
>           DATA(I,J) = DATA(I,J)*MULT
>       enddo
>       enddo
>       ENDIF
>       ELSE
>       read(NF,key=keyA,keyid=0) idat2
> C     TYPE *,(IDAT2.IDAT(70,I),I=1,70,10)
>         do i=1,NLON
>       do j=1,NLAT
>         if (idat2.idat(i,j).gt.-9000.) then 
> C CONVERT AND  MULTIPLY BY FACTOR
>           DATA(I,J) = FLOAT(IDAT2.idat(I,J))*MULT
>       else
>       data(i,j)=-999.
>         endif
>       enddo
>       enddo
>       ENDIF
> c     TYPE *,(DATa(70,I),I=1,70,15)
>       RETURN
>       END
> 
> 
> ------- End of Forwarded Message
>