2016-07-12 34 views
0

我有一個netCDF 4.4文件,它有四個維度 - 時間,緯度,經度,水平。我希望能夠從這個文件中讀取緯度和經度值,並將它們存儲在一維數組中。但是我不知道緯度數組或經度數組的大小,所以我不能在聲明兩個數組lat和lon時分配大小。每個將被處理的netCDF文件將具有它自己的lat和lon大小,所以不可能在編譯時聲明它。所以我用Fortran的分配聲明瞭這一點,但這並沒有達到我想要的效果。它打印0作爲lats和lons的值。 我哪裏錯了?如何在運行時分配Netcdf fortran數組的大小?

Variables for input Z file : height level 
    character*80 in_cfn,varname 
    integer  reason,i,in_ndim,ierr 
    integer ndims_in, nvars_in, ngatts_in, unlimdimid_in 
    integer lat_varid,lon_varid 
    character*(*) LAT_NAME, LON_NAME 
    parameter (LAT_NAME='lat', LON_NAME='lon') 
    real,allocatable, dimension(:) :: lats 
    real,allocatable, dimension(:) :: lons 

    call system('ls hgt_*.nc > hgtFiles.txt') 

    open(10,file='hgtFiles.txt',action="read") 
    varname = "hgt" 
    do 
    read(10,*,IOSTAT=reason) in_cfn 
    if (reason/=0) EXIT 
    print *,in_cfn 
    retval = nf_open(in_cfn,NF_NOWRITE,ncid) 
    if (retval .ne. nf_noerr) call handle_err(retval) 
    retval = nf_inq(ncid,ndims_in,nvars_in,ngatts_in,unlimdimid_in) 
    retval = nf_inq_varid(ncid,LAT_NAME,lat_varid) 
    if (retval .ne. nf_noerr) call handle_err(retval) 
    retval = nf_inq_varid(ncid, LON_NAME, lon_varid) 
    if (retval .ne. nf_noerr) call handle_err(retval) 
    retval = nf_get_var_real(ncid, lat_varid, lats) 
    if (retval .ne. nf_noerr) call handle_err(retval) 
    retval = nf_get_var_real(ncid, lon_varid, lons) 
    print *,size(lons) 

    end do 
    close(10) 

    stop 

    end 

回答

2

我通常使用Fortran 90 NetCDF Module。但Fortran 77 interface工程類似。

您可以用電話獲得尺寸的長度NF_INQ_DIMLEN

我做了一個Fortran 77兼容的(希望)計劃緯度數組中讀取並打印到屏幕上:

 PROGRAM READ_LAT 
     IMPLICIT NONE 
     INCLUDE 'netcdf.inc' 
     INTEGER NCID, LATID, LATVARID, LATLEN 
     REAL LATDATA[ALLOCATABLE](:) 
     CHARACTER*(*) FILENAME 
     CHARACTER*(*) LATNAME 
     PARAMETER (FILENAME='data.nc', LATNAME='lat') 

     CALL CHECK(NF_OPEN(FILENAME, NF_NOWRITE, NCID)) 
     CALL CHECK(NF_INQ_DIMID(NCID, LATNAME, LATID)) 
     CALL CHECK(NF_INQ_DIMLEN(NCID, LATID, LATLEN)) 

     ALLOCATE(LATDATA(LATLEN)) 

     CALL CHECK(NF_INQ_VARID(NCID, LATNAME, LATVARID)) 
     CALL CHECK(NF_GET_VAR_REAL(NCID, LATVARID, LATDATA)) 

     WRITE(*, '(F10.4)') LATDATA 

     CALL CHECK(NF_CLOSE(NCID)) 

     CONTAINS 

     SUBROUTINE CHECK(ERRORCODE) 
      IMPLICIT NONE 
      INTEGER ERRORCODE 
      IF (ERRORCODE .ne. NF_NOERR) THEN 
       PRINT *, "Encountered Error ", ERRORCODE 
       STOP 
      END IF 
     END SUBROUTINE 
     END PROGRAM 

相同的節目的Fortran 90+:

program read_latitude 
    use netcdf 
    implicit none 
    integer :: ncid, latid, latvarid, latlen 
    real, allocatable :: latdata(:) 
    character(len=*), parameter :: filename='data.nc' 
    character(len=*), parameter :: latname = 'lat' 

    call check(nf90_open(filename, NF90_NOWRITE, ncid)) 
    call check(nf90_inq_dimid(ncid, latname, latid)) 
    call check(nf90_inquire_dimension(ncid, latid, len=latlen)) 

    allocate(latdata(latlen)) 

    call check(nf90_inq_varid(ncid, latname, latvarid)) 
    call check(nf90_get_var(ncid, latvarid, latdata)) 

    write(*, '(F10.4)') latdata 

    call check(nf90_close(ncid)) 

    contains 

     subroutine check(errorcode) 
      implicit none 
      integer, intent(in) :: errorcode 
      if (errorcode /= NF90_NOERR) then 
       write(*, '(A, I0)') "Encountered Error ", errorcode 
       write(*, '(A)') nf90_strerror(errorcode) 
       stop 1 
      end if 
     end subroutine check 
end program read_latitude 
+0

感謝您的答覆。假設我確實調用了這個函數,我還需要使用allocate dimension(:)來分配數組嗎? – gansub

+0

是的。這個電話會給你維度的長度。然後你知道你需要分配多大的數組。 – chw21

+0

只是爲了讓你知道,我添加了兩個例子,一個用於F77,另一個用於F90 + – chw21