關於WRF OUTPUT檔的處理

cherryloving

大家好:

小弟最近在使用WRF，WRF OUTPUT檔是NC檔
我寫了一段FORTRAN程式來讀取它的變數，再將變數寫成bin檔

可是發現，不管怎麼讀都讀不到雷達回波及雨資料

想請論壇的高手們幫忙看一下
希望有高手能不吝解答 謝謝

程序碼如附件

readwrf.txt (6.3 KB, 下载次数: 2)

2013-1-28 17:33 上传

点击文件名下载附件

      program ncepdata
      character*80 input

      print*,'Enter input file'
      read(*,'(a)')input
c     input="wrfout_d03_2012-08-12_06:00:00"
      call readwrf(input)

      stop
      end

      subroutine readwrf(input)
c
c  read wrf nc file
c
      parameter (npmax=256*256*65)
      character*80 input,tt*13,out*80
      character*40 var,fname,xtype,name,units,time0,time1
      dimension dat(npmax),dat2(npmax)

      nfopen=nf_open(input,ncnowrit,ncid)
      nfinq=nf_inq(ncid,ndim,nvar,ngatte,ire)
c     print*,'ndim=',ndim,' nvar=',nvar

      do i=1,ndim
         name=' '
         inqdim=nf_inq_dim(ncid,i,name,len)
         if(name.eq.'west_east')nx=len
         if(name.eq.'south_north')ny=len
         if(name.eq.'bottom_top')nz=len
         if(name.eq.'west_east_stag')nx1=len
         if(name.eq.'south_north_stag')ny1=len
         if(name(1:15).eq.'bottom_top_stag')nz1=len
         print*,i,' dim=',name(1:20),len
      enddo
      print*,nx,ny,nz,nx1,ny1,nz1

      do i=1,nvar
         fname=' '
         units=' '
         name=' '
         inqvar=nf_inq_var(ncid,i,fname,xtype,idim,idimid,iatt)
         inqdim=nf_inq_dimlen(ncid,i,ilen)
         igtatt=nf_get_att_text(ncid,i,'units',units)
         igtatt=nf_get_att_text(ncid,i,'long_name',name)
c        print*,'var=',fname(1:10),idim,'  ',name(1:20),units(1:26)
      enddo

      time0=' '
      time1=' '
      var='Times'
      inqvar=nf_inq_varid(ncid,var,idvar)
      igtvar=nf_get_var_text(ncid,idvar,time1)

      igtatt=nf_get_att_text(ncid,0,'SIMULATION_START_DATE',time0)
      igtatt=nf_get_att_real(ncid,0,'DX',dx)
      igtatt=nf_get_att_real(ncid,0,'DY',dy)
      igtatt=nf_get_att_real(ncid,0,'DT',dt)
      igtatt=nf_get_att_real(ncid,0,'CEN_LAT',rlat0)
      igtatt=nf_get_att_real(ncid,0,'CEN_LON',rlon0)
      igtatt=nf_get_att_real(ncid,0,'TRUELAT1',tlat1)
      igtatt=nf_get_att_real(ncid,0,'TRUELAT2',tlat2)

      read(time0,11)iy0,im0,id0,ih0,in0,is0
      read(time1,11)iy1,im1,id1,ih1,in1,is1
  11  format(i4,5(1x,i2))
      write(tt,'(i4.4,2i2.2,a1,2i2.2)')iy1,im1,id1,'_',ih1,in1

      nz1=nz+1
      print*,'t0=',time0
      print*,'t1=',time1,' dt(min)=',dt
      print*,'grid=',nx,ny,nz,nx1,ny1,nz1,dx,dy
      print*,'map=',rlat0,rlon0,tlat1,tlat2

      var='ZNU'
      call getdata(ncid,var,dat,nz,1,1)
      var='ZNW'
      call getdata(ncid,var,dat,nz1,1,1)

      var='U'
      out=tt//'.u.bin'
      call getdata(ncid,var,dat,nx1,ny,nz)
      call udat(dat,nx1,ny,nz,dat2,nx)
      call output(dat2,nx,ny,nz,out)

      var='V'
      out=tt//'.v.bin'
      call getdata(ncid,var,dat,nx,ny1,nz)
      call vdat(dat,nx,ny1,nz,dat2,ny)
      call output(dat2,nx,ny,nz,out)

      var='W'
      out=tt//'.w.bin'
      call getdata(ncid,var,dat,nx,ny,nz1)
      call wdat(dat,nx,ny,nz1,dat2,nz)
      call output(dat2,nx,ny,nz,out)

      var='PH'
      out=tt//'.z.bin'
      call getdata(ncid,var,dat,nx,ny,nz1)
      var='PHB'
      call getdata(ncid,var,dat2,nx,ny,nz1)
      call add(dat,dat2,nx,ny,nz1)
      call wdat(dat,nx,ny,nz1,dat2,nz)
      call output(dat,nx,ny,nz,out)

      var='T'
      call getdata(ncid,var,dat,nx,ny,nz)

      var='P'
      out=tt//'.p.bin'
      call getdata(ncid,var,dat,nx,ny,nz)
      var='PB'
      call getdata(ncid,var,dat2,nx,ny,nz)
      call add(dat,dat2,nx,ny,nz)
      call output(dat,nx,ny,nz,out)

      var='T2'
      out=tt//'.t2.bin'
      call getdata(ncid,var,dat,nx,ny,1)
      call output(dat,nx,ny,1,out)

      var='PSFC'
      call getdata(ncid,var,dat,nx,ny,1)

      var='U10'
      out=tt//'.u10.bin'
      call getdata(ncid,var,dat,nx,ny,1)
      call output(dat,nx,ny,1,out)

      var='V10'
      out=tt//'.v10.bin'
      call getdata(ncid,var,dat,nx,ny,1)
      call output(dat,nx,ny,1,out)

      var='HGT' ! terrain height
      call getdata(ncid,var,dat,nx,ny,1)

      var='RAINC'
      out=tt//'.rain.bin'
      call getdata(ncid,var,dat,nx,ny,1)
      var='RAINSH'
      call getdata(ncid,var,dat2,nx,ny,1)
      call add(dat,dat2,nx,ny,1)
      var='RAINNC'
      call getdata(ncid,var,dat2,nx,ny,1)
      call add(dat,dat2,nx,ny,1)
      call output(dat,nx,ny,1,out)

      var='REFL_10CM'
      out=tt//'.dbz.bin'
      call getdata(ncid,var,dat,nx,ny,nz)
      call output(dat,nx,ny,nz,out)

      nfclose=nf_close(ncid)

      return
      end

      subroutine getdata(ncid,var,dat,nx,ny,nz)
c
c  get variable data
c
      character*31 var
      dimension istart(10),icount(10),dat(nx*ny*nz)

      istart(1)=1
      istart(2)=1
      istart(3)=1
      istart(4)=1
      icount(1)=nx
      icount(2)=ny
      icount(3)=nz
      icount(4)=1
      inqvar=nf_inq_varid(ncid,var,idvar)
      igtvar=nf_get_vara_real(ncid,idvar,istart,icount,dat)

      npt=nx*ny*nz
      dmin=1.e20
      dmax=-dmin
      do i=1,npt
         dmin=min(dmin,dat(i))
         dmax=max(dmax,dat(i))
      enddo
      print*,nx,ny,nz,' ',var(1:20),dmin,dmax

      subroutine udat(dat,nx1,ny,nz,dat2,nx)
      dimension dat(nx1,ny,nz),dat2(nx,ny,nz)
      do k=1,nz
      do j=1,ny
      do i=1,nx
         dat2(i,j,k)=(dat(i,j,k)+dat(i+1,j,k))/2.
      enddo
      enddo
      enddo
      return
      end

      subroutine vdat(dat,nx,ny1,nz,dat2,ny)
      dimension dat(nx,ny1,nz),dat2(nx,ny,nz)
      do k=1,nz
      do j=1,ny
      do i=1,nx
         dat2(i,j,k)=(dat(i,j,k)+dat(i,j+1,k))/2.
      enddo
      enddo
      enddo
      return
      end

      subroutine wdat(dat,nx,ny,nz1,dat2,nz)
      dimension dat(nx,ny,nz1),dat2(nx,ny,nz)
      do k=1,nz
      do j=1,ny
      do i=1,nx
         dat2(i,j,k)=(dat(i,j,k)+dat(i,j,k+1))/2.
      enddo
      enddo
      enddo
      return
      end

      subroutine output(dat,nx,ny,nz,out)
      dimension dat(nx*ny*nz)
      character*80 out
      np=nx*ny*nz
      open(51,file=out,form='unformatted',status='unknown',
     +        access='direct',recl=np*4)
      write(51,rec=1)(dat(i),i=1,np)
      close(51)
      return
      end


      subroutine add(dat,dat2,nx,ny,nz)
      dimension dat(nx,ny,nz),dat2(nx,ny,nz)
      do k=1,nz
      do j=1,ny
      do i=1,nx
         dat(i,j,k)=dat(i,j,k)+dat2(i,j,k)
      enddo
      enddo
      enddo
      return
      end

三色堇

我也遇到了一些问题，用nf_open始终没法打开wrfout的文件