module initInterp2dvar ! ------------------------------------------------------------------------ ! !DESCRIPTION: ! This module defines a class that holds information about a 2-d variable, to ! facilitate reading & writing ! ! !USES: #include "shr_assert.h" use abortutils , only: endrun use shr_log_mod , only: errMsg => shr_log_errMsg use shr_kind_mod , only: r8 => shr_kind_r8, i4=>shr_kind_i4 use clm_varctl , only: iulog use ncdio_pio , only: file_desc_t, var_desc_t, ncd_inqvid, ncd_inqvdids use ncdio_pio , only: ncd_inqdlen, ncd_inqdname, ncd_io, ncd_getatt use pio , only: pio_inq_varid, pio_get_var implicit none private save ! Public types public :: interp_2dvar_type type :: interp_2dvar_type private character(len=:), allocatable :: varname ! variable name type(var_desc_t) :: vardesc ! variable descriptor on file integer :: varid ! variable ID on file type(file_desc_t), pointer :: ncid ! pointer to netcdf ID corresponding to this variable logical :: file_is_dest ! true if this is on the dest file, false if on the source file character(len=16) :: vec_dimname ! dimension name of vector dimension (e.g., 'col') character(len=16) :: lev_dimname ! dimension name of level dimension integer :: vec_beg ! beginning index of vector dimension integer :: vec_end ! ending index of vector dimension integer :: nlev ! size of level dimension logical :: switchdim ! true if dimensions are 'switched' for this variable contains ! Public routines generic :: readvar => readvar_int, readvar_double !TYPE int,double procedure :: readvar_{TYPE} generic :: writevar => writevar_int, writevar_double !TYPE int,double procedure :: writevar_{TYPE} generic :: readlevel => readlevel_int, readlevel_double ! read one level !TYPE int,double procedure :: readlevel_{TYPE} procedure :: get_varname ! get variable name procedure :: get_vec_dimname ! get dimension name of vector dimension (e.g., 'col') procedure :: get_lev_dimname ! get dimension name of level dimension procedure :: get_vec_beg ! get beginning index of vector dimension procedure :: get_vec_end ! get ending index of vector dimension procedure :: get_vec_npts ! get number of points in vector dimension procedure :: get_nlev ! get number of levels ! Private routines procedure, private :: set_switchdim end type interp_2dvar_type interface interp_2dvar_type module procedure constructor end interface interp_2dvar_type character(len=*), parameter, private :: sourcefile = & __FILE__ contains ! ======================================================================== ! Constructors ! ======================================================================== !----------------------------------------------------------------------- function constructor(varname, ncid, file_is_dest, bounds) & result(this) ! ! !DESCRIPTION: ! Creates an interp_2dvar_type object ! ! !USES: use initInterpBounds, only : interp_bounds_type ! ! !ARGUMENTS: type(interp_2dvar_type) :: this ! function result character(len=*) , intent(in) :: varname ! variable name type(file_desc_t) , target, intent(inout) :: ncid ! netcdf id logical , intent(in) :: file_is_dest ! true if this is the dest file, false if the source file type(interp_bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: integer :: vec_dimnum ! dimension number of the vector dimension (e.g., 'col') integer :: lev_dimnum ! dimension number of the level dimension integer :: dimids(2) ! dimension IDs character(len=*), parameter :: subname = 'interp_2dvar_type constructor' !----------------------------------------------------------------------- this%varname = varname this%ncid => ncid this%file_is_dest = file_is_dest call ncd_inqvid(this%ncid, trim(varname), this%varid, this%vardesc) call this%set_switchdim() ! ------------------------------------------------------------------------ ! Get information related to dimensions ! ------------------------------------------------------------------------ if (this%switchdim) then lev_dimnum = 1 vec_dimnum = 2 else lev_dimnum = 2 vec_dimnum = 1 end if call ncd_inqvdids(this%ncid, dimids, this%vardesc) ! Get size of level dimension call ncd_inqdlen(this%ncid, dimids(lev_dimnum), this%nlev) ! Get name of level dimension call ncd_inqdname(this%ncid, dimids(lev_dimnum), this%lev_dimname) ! Get name of vector dimension call ncd_inqdname(this%ncid, dimids(vec_dimnum), this%vec_dimname) this%vec_beg = bounds%get_beg(this%vec_dimname) this%vec_end = bounds%get_end(this%vec_dimname) end function constructor ! ======================================================================== ! Public methods ! ======================================================================== !----------------------------------------------------------------------- !TYPE int,double subroutine readvar_{TYPE}(this, data) ! ! !DESCRIPTION: ! Reads variable from file. ! ! The 'data' variable is allocated here, and must be deallocated by the caller. ! ! !USES: ! ! !ARGUMENTS: ! 'this' is intent(inout) for the sake of ncid class(interp_2dvar_type), intent(inout) :: this {VTYPE}, pointer, intent(out) :: data(:,:) ! [vec, lev] ! ! !LOCAL VARIABLES: {VTYPE}, pointer :: data_transpose(:,:) ! transpose of data array if switchdim is true character(len=*), parameter :: subname = 'readvar' !----------------------------------------------------------------------- if (.not. this%file_is_dest) then ! NOTE(wjs, 2015-11-22) The code is set up to be able to read properly from the ! source file, and this has been tested. However, for relatively high resolution ! source files (f09 or higher, and/or higher-than-standard number of vertical ! levels) it can cause out-of-memory problems on yellowstone that bring the whole ! system to a temporary stand-still. So, to be safe, we currently abort if this is ! attempted. write(iulog,*) subname//' ERROR: Attempt to read all levels at once from source file.' write(iulog,*) 'Although this is possible in theory, it is generally a bad idea' write(iulog,*) 'because it can result in extreme memory usage.' write(iulog,*) 'In some cases, this can cause system hangs that affect all users.' write(iulog,*) 'If you want to try anyway, you can remove this endrun call.' call endrun(msg='Attempt to read all levels at once from source file '// & errMsg(sourcefile, __LINE__)) end if allocate(data(this%vec_beg:this%vec_end, this%nlev)) ! The destination file uses the model's decomposition, whereas the source file does ! not. Thus, we need to handle the reads differently. if (this%file_is_dest) then call ncd_io(ncid=this%ncid, varname=trim(this%varname), flag='read', data=data, & dim1name=trim(this%vec_dimname), switchdim=this%switchdim) else if (this%switchdim) then allocate(data_transpose(this%nlev, this%vec_beg:this%vec_end)) call ncd_io(ncid=this%ncid, varname=trim(this%varname), flag='read', & data=data_transpose) data = transpose(data_transpose) deallocate(data_transpose) else call ncd_io(ncid=this%ncid, varname=trim(this%varname), flag='read', & data=data) end if end if end subroutine readvar_{TYPE} !----------------------------------------------------------------------- !TYPE int,double subroutine writevar_{TYPE}(this, data) ! ! !DESCRIPTION: ! Writes variable to file. ! ! Currently only works for the destination file, not the source file. ! ! !USES: ! ! !ARGUMENTS: ! 'this' is intent(inout) for the sake of ncid class(interp_2dvar_type), intent(inout) :: this {VTYPE}, pointer, intent(in) :: data(:,:) ! [vec, lev] ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'writevar' !----------------------------------------------------------------------- if (this%file_is_dest) then call ncd_io(ncid=this%ncid, varname=trim(this%varname), flag='write', data=data, & dim1name=trim(this%vec_dimname), switchdim=this%switchdim) else ! NOTE(wjs, 2015-10-13) In principle, we could probably handle writes to the ! source file the same way as we handle reads. But since this is untested, I'm ! currently aborting in this case. call endrun(msg='ERROR: unhandled attempt to write variable to source file '// & errMsg(sourcefile, __LINE__)) end if end subroutine writevar_{TYPE} !----------------------------------------------------------------------- !TYPE int,double subroutine readlevel_{TYPE}(this, data, level) ! ! !DESCRIPTION: ! Reads one level of variable from file. ! ! In contrast to readvar, this routine does NOT allocate the data: it expects data to ! already be allocated of the correct size. (This is to avoid requiring repeated ! allocations and deallocations when reading one level at a time.) ! ! !USES: ! ! !ARGUMENTS: ! 'this' is intent(inout) for the sake of ncid class(interp_2dvar_type), intent(inout) :: this {VTYPE}, pointer :: data(:) integer, intent(in) :: level ! level index to read ! ! !LOCAL VARIABLES: integer :: status ! netCDF return code type(Var_desc_t) :: vardesc ! pio variable descriptor integer :: starts(2) integer :: counts(2) character(len=*), parameter :: subname = 'readlevel' !----------------------------------------------------------------------- SHR_ASSERT_ALL((lbound(data) == (/this%vec_beg/)), errMsg(sourcefile, __LINE__)) SHR_ASSERT_ALL((ubound(data) == (/this%vec_end/)), errMsg(sourcefile, __LINE__)) if (this%file_is_dest) then ! NOTE(wjs, 2015-11-22) As far as I can tell, ncdio_pio currently doesn't provide ! a robust mechanism for reading in a single level. Since we don't actually need ! this functionality right now, I am aborting rather than trying to either (a) add ! the functionality to ncdio_pio, or (b) handle it here with direct calls to pio. call endrun(msg='ERROR: unhandled attempt to read one level from dest file '//& errMsg(sourcefile, __LINE__)) else ! TODO(wjs, 2015-11-22) I don't like that we're making direct pio calls here. ! Ideally we'd go through ncdio_pio, but currently that module doesn't provide the ! necessary functionality. status = pio_inq_varid(this%ncid, trim(this%varname), vardesc) if (this%switchdim) then starts(1) = level counts(1) = 1 starts(2) = 1 counts(2) = this%get_vec_npts() else starts(1) = 1 counts(1) = this%get_vec_npts() starts(2) = level counts(2) = 1 end if status = pio_get_var(this%ncid, vardesc, starts, counts, data) end if end subroutine readlevel_{TYPE} function get_varname(this) result(varname) ! Get variable name character(len=:), allocatable :: varname ! function result class(interp_2dvar_type), intent(in) :: this varname = this%varname end function get_varname function get_vec_dimname(this) result(vec_dimname) ! Get name of vector dimension (e.g., 'col') character(len=:), allocatable :: vec_dimname ! function result class(interp_2dvar_type), intent(in) :: this vec_dimname = this%vec_dimname end function get_vec_dimname function get_lev_dimname(this) result(lev_dimname) ! Get name of level dimension character(len=:), allocatable :: lev_dimname ! function result class(interp_2dvar_type), intent(in) :: this lev_dimname = this%lev_dimname end function get_lev_dimname integer function get_vec_beg(this) ! Get beginning index of vector dimension class(interp_2dvar_type), intent(in) :: this get_vec_beg = this%vec_beg end function get_vec_beg integer function get_vec_end(this) ! Get ending index of vector dimension class(interp_2dvar_type), intent(in) :: this get_vec_end = this%vec_end end function get_vec_end integer function get_vec_npts(this) ! Get number of points in vector dimension class(interp_2dvar_type), intent(in) :: this get_vec_npts = (this%vec_end - this%vec_beg + 1) end function get_vec_npts integer function get_nlev(this) ! Get number of levels class(interp_2dvar_type), intent(in) :: this get_nlev = this%nlev end function get_nlev ! ======================================================================== ! Private methods ! ======================================================================== !----------------------------------------------------------------------- subroutine set_switchdim(this) ! ! !DESCRIPTION: ! Sets the switchdim attribute by reading the appropriate metadata from the file ! ! !USES: use restUtilMod , only: iflag_noswitchdim, iflag_switchdim ! ! !ARGUMENTS: class(interp_2dvar_type), intent(inout) :: this ! ! !LOCAL VARIABLES: integer :: switchdim_flag character(len=*), parameter :: subname = 'set_switchdim' !----------------------------------------------------------------------- call ncd_getatt(this%ncid, this%varid, 'switchdim_flag', switchdim_flag) select case (switchdim_flag) case (iflag_switchdim) this%switchdim = .true. case (iflag_noswitchdim) this%switchdim = .false. case default write(iulog,*) subname//' ERROR: unknown switchdim_flag: ', switchdim_flag call endrun('ERROR: unknown switchdim flag '//errMsg(sourcefile, __LINE__)) end select end subroutine set_switchdim end module initInterp2dvar