FORTRAN程序读取ctl文件,并写成netcdf4格式

wsdwsd

!读取grads的ctl描述的数据,生成netcdf4格式
!编译方法: 请根据自己实际环境进行修改
!pgf90  -I/opt/local/netcdf/include   -L/opt/local/netcdf/lib -lnetcdff  nc4write3D.f90  -o nc4write3D.exe
!by wsd@cma.gov.cn

      !nc4相关定义
      use netcdf
      parameter (NDIMS2=2,NDIMS3=3,maxitem=20)
      integer iddimsize2(NDIMS2),iddimsize3(NDIMS3),chunks(NDIMS2),chunks3(NDIMS3),idnc,id_lat,id_lon,id_lev,iddim_lat,iddim_lon,idallvar(maxitem)
      real, dimension(:), allocatable :: lat_array,lon_array,lev_array
      integer start(NDIMS3), count(NDIMS3)

      !数据相关的定义
      integer ,iy0,ix0,iz0,ctlvarflag(maxitem)
      real, dimension(:,:), allocatable :: realvar
      real, dimension(:,:,:), allocatable :: realvar3d
      real levels(1000)
      character*20 chvarname(2,maxitem)  !1:varname,2:unitname
      character*180 cmd,filectl,nameofdat,nameout

      numarg = iargc()
      call getarg(0,cmd)
      if(numarg<1)then  ! .or. numarg /=2
         write(*,'(a)')'Usage: '//trim(cmd)//' xxxx.ctl'
         write(*,'(a)')'   ie: '//trim(cmd)//' shelt.ctl'
         write(*,'(a)')'   ie: '//trim(cmd)//' latlon_shelt.ctl'
         stop
      endif
      call getarg(1,filectl)
      call readctl(filectl,nameofdat,nameout,ctlvarflag,ix0,iy0,iz0,startx,starty,stepx,stepy,nvar,chvarname,maxitem,undef,if3d,levels,1000)   
      ifgeog=1;if(abs(startx-starty-stepx+stepy)<0.01)ifgeog=0  !检测ctl是否有地理信息
      !write(*,*)abs(startx-starty-stepx+stepy),'====',ifgeog
      ALLOCATE(lon_array(ix0),lat_array(iy0),lev_array(iz0),realvar(ix0,iy0))
      if(if3d==1)then
        ALLOCATE(realvar3d(ix0,iy0,iz0))
        !write(*,*)isz0,' levels:',levels(1:iz0)
      endif
      open(10,file=nameofdat,status='old',form='binary',err=999)
      iret = nf90_create(trim(nameout)//'.nc', nf90_netcdf4, idnc)  ! Create the file.
      
      !设置各坐标变量
      if(ifgeog==1)then!是经纬度格式,则设置地图经纬度相关信息
        write(*,'(4i5,4f9.3,f10.1,i2)')ix0,iy0,iz0,nvar,startx,starty,stepx,stepy,undef,if3d
        do i=1,ix0
           lon_array(i)=startx+real(i-1)*stepx
        enddo
        do j=1,iy0
          lat_array(j)=starty+real(j-1)*stepy
        enddo
        do j=1,iz0
          lev_array(j)=levels(j)
        enddo
        !定义维数
        iret = nf90_def_dim(idnc, 'latitude',  iy0, iddim_lat);if (iret .ne. nf90_noerr) call handle_err(iret)
        iret = nf90_def_dim(idnc, 'longitude', ix0, iddim_lon);if (iret .ne. nf90_noerr) call handle_err(iret)
        iret = nf90_def_dim(idnc, 'level',     iz0, iddim_lev);if (iret .ne. nf90_noerr) call handle_err(iret)
        !定义经纬度变量
        iret = nf90_def_var(idnc, 'latitude',  NF90_REAL, iddim_lat,id_lat);if (iret .ne. nf90_noerr) call handle_err(iret)
        iret = nf90_def_var(idnc, 'longitude', NF90_REAL, iddim_lon,id_lon);if (iret .ne. nf90_noerr) call handle_err(iret)
        iret = nf90_def_var(idnc, 'level',     NF90_REAL, iddim_lev,id_lev);if (iret .ne. nf90_noerr) call handle_err(iret)
        !写经纬度属性
        iret = nf90_put_att(idnc, id_lat, 'units', 'degrees_north');if (iret .ne. nf90_noerr) call handle_err(iret)
        iret = nf90_put_att(idnc, id_lon, 'units', 'degrees_east');if (iret .ne. nf90_noerr) call handle_err(iret)
        if(if3d==1)then
          iret = nf90_put_att(idnc, id_lev, 'units', 'm');      if (iret .ne. nf90_noerr) call handle_err(iret)
        endif
      else !无地理属性,直接写维数
        write(*,'(4i5,f10.1)')ix0,iy0,iz0,nvar,undef
        iret = nf90_def_dim(idnc, 'x', ix0, iddim_lon)
        iret = nf90_def_dim(idnc, 'y', iy0, iddim_lat)
        iret = nf90_def_dim(idnc, 'z', iz0, iddim_lev)
      endif
      iddimsize2 = (/ iddim_lon, iddim_lat /)
      iddimsize3 = (/ iddim_lon, iddim_lat, iddim_lev /)
      chunks = (/ix0,iy0/)
      chunks3 = (/ix0,iy0,1/)  !对于大的3D数组,只能这样设置,较小的可以将Z方向也设上
      levelofdeflate=1         !压缩强度,1已经达到足够的压缩比了,9则耗费太多的时间

      !设置各要素变量
      do j=1,nvar  !依次定义各变量
        write(*,'(a,i2,2x,a)')'defing ',j,trim(chvarname(1,j))//'  '//trim(chvarname(2,j))
        if(ctlvarflag(j)==2)then
          iret = nf90_def_var(idnc, trim(chvarname(1,j)),NF90_REAL,iddimsize2,idallvar(j),chunksizes=chunks)!,shuffle=.TRUE.,deflate_level=deflate_level)        
        elseif(ctlvarflag(j)==3)then
          iret = nf90_def_var(idnc, trim(chvarname(1,j)),NF90_REAL,iddimsize3,idallvar(j),chunksizes=chunks3)!,shuffle=.TRUE.,deflate_level=deflate_level)        
        endif
        iret = nf90_def_var_deflate(idnc, idallvar(j), 1, 1, levelofdeflate)
        iret = nf90_put_att(idnc, idallvar(j), 'units',trim(chvarname(2,j)));if (iret .ne. nf90_noerr) call handle_err(iret)  ! len_trim(chvarname(1,j)),
        iret = nf90_put_att(idnc, idallvar(j), 'missing_value', undef)
      enddo
      iret = nf90_enddef(idnc);if (iret .ne. nf90_noerr) call handle_err(retval)  !完成定义工作

      !写入各数据
      if(ifgeog==1)then!是经纬度格式,非经纬度的3D数组则不用设置
        iret = nf90_put_var(idnc, id_lat, lat_array)
        iret = nf90_put_var(idnc, id_lon, lon_array)
        if(if3d==1)iret = nf90_put_var(idnc, id_lev, lev_array)
      endif

      count=(/ix0,iy0,iz0/);start=(/1,1,1/)  !3D和4D数据写入必须要这些
      do j=1,nvar
        if(ctlvarflag(j)==2)then
          write(*,*)'writing 2D',j
          read(10)realvar
          iret = nf90_put_var(idnc, idallvar(j), realvar);if (iret .ne. nf90_noerr) call handle_err(iret) !写入各变量
        elseif(ctlvarflag(j)==3)then
          read(10)realvar3d
          write(*,*)'writing 3D',j
          iret = nf90_put_var(idnc, idallvar(j), realvar3d,start=start,count=count);if (iret .ne. nf90_noerr) call handle_err(iret) !写入各变量
        endif
      enddo
      close(10)
      iret = nf90_close(idnc);if (iret .ne. nf90_noerr) call handle_err(iret)
      write(*,*)'SUCCESS writing to NC file: '//trim(nameout)//'.nc'
      stop ''
  999 write(*,*)'Not found dat file:'//trim(nameofdat)//'=='
      end

      subroutine handle_err(errcode)
      use netcdf
      implicit none
      !include 'netcdf.inc'
      integer errcode
      print *, 'Error: ', nf90_strerror(errcode)
      stop 2
      end

      subroutine readctl(filectl,nameofdat,nameout,ctlvarflag,ix0,iy0,iz0,startx,starty,stepx,stepy,ivar,chvarname,maxitem,undef,if3d,levels,maxlev)   
      character*(*) filectl,nameofdat,nameout,chvarname(2,maxitem),ch*200,ch1*30
      real levels(maxlev)
      integer ctlvarflag(maxitem)
      if3d=0  !作为是否有3D数据的标志
      ctlvarflag=-1
      open(100,file=filectl,status='old')
      read(100,'(a200)')ch
      ilen=len_trim(ch)
      if(ch(ilen:ilen)<'!')ch(ilen:ilen)=''  !避免windows格式带来的回车换行带来的影响        
      nameofdat=ch(6:len_trim(ch))  !'dset ./AN...' 防止有/误作转义符
      if(nameofdat(1:1)=='^')nameofdat=nameofdat(2:len_trim(nameofdat)) !忽略grads里的特有符号
  100 read(100,'(a200)')ch             !
      if(ch(1:5)=='title' .or. ch(1:5)=='TITLE')goto 100     !忽略title
      read(ch,*)ch1,undef     
      read(100,*)ch,ix0,ch,startx,stepx     
      read(100,*)ch,iy0,ch,starty,stepy     
      read(100,*)ch,iz0,ch,(levels(i),i=1,iz0)  !可能按行排列,或本行随后排放
      read(100,*)                               !tdef  暂时忽略时间序列
      read(100,*)ch,ivar
      if(ivar>maxitem)ivar=maxitem
      do i=1,ivar
        read(100,*)chvarname(1,i),i1,i2,chvarname(2,i)
        if(i1/=0)then !三维场
          ctlvarflag(i)=3
          if3d=1
        else
          ctlvarflag(i)=2  !二维场
        endif
      enddo   
      close(100)
      nameout=trim(nameofdat)  !准备输出文件名,抛弃原路径和后缀
      i=index(nameout,'/',BACK=.TRUE.)  !抛弃原路径
      if(i>0)nameout=nameout(i+1:200)
      !write(*,*)nameout(i+1:i+10)
      j=index(nameout,'.',BACK=.TRUE.)  !抛弃后缀
      if(j>0)nameout(j:200)=''
      !write(*,*)ivar,'=='//trim(nameofdat)//'==>'//trim(nameout)//'=='
      do i=1,ivar
        write(*,'(2i4,a,a20,2x,a)')i,ctlvarflag(i),'D',trim(chvarname(1,i)),trim(chvarname(2,i))
      enddo   
      end      

又是那隻貓

谢谢分享

ningking1986

感谢分享啊~！

sunflower

感谢楼主分，楼主好人啊，还不需要花钱

cibaouc

楼主，netcdf4的读取怎么读取的