2015-01-01 24 views
0

我试图修改一个Fortran 90代码,它以NetCDF经典格式将二维数组写入输出。我希望变量有一个额外的时间维度(即,它将是一个3D变量),在模型的积分时间内的每个相应的时间步长上打印它。在NetCDF中写入作为时间函数的变量

我不确定它是如何完成的;我很欣赏任何建议,尽可能有效地做到这一点(也是最小的文件大小)。

subroutine writenetcdffile(array,argtitle) 
    use netcdf 
    implicit none 
    real, intent(IN), dimension(:,:) :: array 
    character*(*),intent(IN) :: argtitle 

    integer :: file_id, xdim_id, ydim_id 
    integer :: array_id 
    integer, dimension(2) :: arrdims 
! character(len=*) :: argtitle = Flag_in 

    integer :: i, j 
    integer :: ierr 

    i = size(array,1) 
    j = size(array,2) 

    ! create the file 
    ierr = nf90_create(path='test.nc', cmode=NF90_CLOBBER, ncid=file_id) 

    ! define the dimensions 
    ierr = nf90_def_dim(file_id, 'X', i, xdim_id) 
    ierr = nf90_def_dim(file_id, 'Y', j, ydim_id) 

    ! now that the dimensions are defined, we can define variables on them,... 
    arrdims = (/ xdim_id, ydim_id /) 
    ierr = nf90_def_var(file_id, 'Array', NF90_REAL, arrdims, array_id) 

    ! ...and assign units to them as an attribute 
    ierr = nf90_put_att(file_id, array_id, "title", argtitle) 

    ! done defining 
    ierr = nf90_enddef(file_id) 

    ! Write out the values 
    ierr = nf90_put_var(file_id, array_id, array) 

    ! close; done 
    ierr = nf90_close(file_id) 
    return 
    end subroutine writenetcdffile 

MODULE Module_NetCDF 

use netcdf 
IMPLICIT NONE 

integer :: file_id, xdim_id, ydim_id, tdim_id 
integer :: array_id(5) 
integer, dimension(3) :: arrdims 

integer :: i, j 
integer :: ierr 

CONTAINS 

SUBROUTINE NetCDF_Init(ICase) 

    IMPLICIT NONE 

    INTEGER :: ICase 

    SELECT CASE(ICase) 
    Case(1) 
     ! create the file 
     ierr = nf90_create(path='test.nc', cmode = NF90_CLOBBER, ncid = file_id) 
    Case(2) 
     ! Reopen the file for writing 
     ierr = nf90_open(path = "test.nc", mode = nf90_write, ncid = file_id) 
     if (ierr /= nf90_noerr) call check(ierr) 
    Case(3) 
     ! close; done 
     ierr = nf90_close(file_id) 
    END SELECT 

    RETURN 
    END SUBROUTINE NetCDF_Init 


    SUBROUTINE NetCDF_Def(Array,ArrayTitle,ArrayUnits) 

    IMPLICIT NONE 

    real, intent(IN), dimension(:,:) :: Array 
    character(*),intent(IN) :: ArrayTitle(5) 
    character(*),intent(IN) :: ArrayUnits(5) 

! Locals 
    integer :: k 

    i = size(Array,1) 
    j = size(Array,2) 

! CALL NetCDF_Init(1) 

    ! define the dimensions 
    ierr = nf90_def_dim(file_id, 'X', i, xdim_id) 
    ierr = nf90_def_dim(file_id, 'Y', j, ydim_id) 
    ierr = nf90_def_dim(file_id, 'Time', nf90_unlimited, tdim_id) 

    ! now that the dimensions are defined, we can define variables on them,... 
    arrdims = (/ xdim_id, ydim_id, tdim_id /) 
    do k = 1,size(ArrayTitle) 
     ierr = nf90_def_var(file_id, ArrayTitle(k), NF90_REAL, arrdims, array_id(k)) 

    ! ...and assign units to them as an attribute 
     ierr = nf90_put_att(file_id, array_id(k), "Units", ArrayUnits(k)) 
    enddo 

    ! done defining 
    ierr = nf90_enddef(file_id) 

    RETURN 
    END SUBROUTINE NetCDF_Def 
SUBROUTINE NetCDF_Write(Array,FlagTitle,NTime) 

    IMPLICIT NONE 

    real, intent(IN), dimension(:,:) :: Array 
    integer,intent(IN) :: NTime 
    character(*),intent(in) :: FlagTitle 

! Locals 
    integer :: J_id 

    IF(FlagTitle.EQ.'ONECOND')THEN 
     J_id = 1 
    ELSEIF(FlagTitle.EQ.'MELTING')THEN 
     J_id = 2 
    ELSEIF(FlagTitle.EQ.'FREEZ_NEW')THEN 
     J_id = 3 
    ELSEIF(FlagTitle.EQ.'TFREEZ')THEN 
     J_id = 4 
    ELSEIF(FlagTitle.EQ.'DFREEZ')THEN 
     J_id = 5  
    ENDIF 

    CALL NetCDF_Init(2) 

    ierr = nf90_put_var(file_id, array_id(j_id), Array, start=[1,1,ntime], count=[i,j,1]) 

    CALL NetCDF_Init(3) 


    RETURN 
    END SUBROUTINE 
SUBROUTINE check(status) 

    IMPLICIT NONE 
    integer, intent (in) :: status 

    IF(status /= nf90_noerr) THEN 
     PRINT *, trim(nf90_strerror(status)) 
     STOP 2 
    ENDIF 
    END SUBROUTINE check 

END MODULE Module_NetCDF 
+1

你给的例程是2D数组......你尝试过3D阵列的是什么?你得到什么错误?就像现在一样,你的问题太广泛了! –

回答

1

你需要做的是确定nf90_unlimited长的时间维度。这将允许您一次将一个二维数组写入一个三维数组,并使该数组的长度未指定。使用startcount可选虚拟参数调用nf90_put_var来指定写入2-d片段的位置。

! create the file 
ierr = nf90_create(path='test.nc', cmode=NF90_CLOBBER, ncid=file_id) 

! define the dimensions 
ierr = nf90_def_dim(file_id, 'X', i, xdim_id) 
ierr = nf90_def_dim(file_id, 'Y', j, ydim_id) 
ierr = nf90_def_dim(file_id, 'Time', nf90_unlimited, tdim_id) 

! now that the dimensions are defined, we can define variables on them,... 
arrdims = (/ xdim_id, ydim_id, tdim_id /) 
ierr = nf90_def_var(file_id, 'Array', NF90_REAL, arrdims, array_id) 

! done defining 
ierr = nf90_enddef(file_id) 

! Time loop 
do n = 1,nm 

    ! Calculations go here 

    ! Write out the values  
    ierr = nf90_put_var(file_id, array_id, array, start=[1,1,n], count=[i,j,1]) 

enddo 

我在我的大多数程序做的是创建该文件,并在一开始定义维度和变量,后来写在一个循环中的字段。如果您的模拟需要很长时间,并且您希望能够在模拟过程中查看输出,请在模型求解器的do-loop中执行打开/写入/关闭步骤。

+0

谢谢@milancurcic。这正是我所需要的,我会试一试。在NetCDFIf中,如果我想在运行时观察输出,我必须在用MATLAB分析之前关闭文件吗? – Jacob

+0

顺便说一句 - 哪一个Call命令对应于每次'nf90_put_var'调用时打开文件? (我假设我不应该使用'nf90_create'重复)。 – Jacob

+0

@Jacob分别使用'nf90_open'和'nf90_close'打开和关闭文件。我相信你需要关闭该文件才能阅读它,但是,我注意到ncview可以在NetCDF文件处于打开状态并正在写入时读取这些文件。试一试。另外,请在线查看Fortran NetCDF用户指南。这非常有用。 – milancurcic