Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

diag_axis io #982

Merged
merged 6 commits into from
Jun 3, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions diag_manager/diag_data.F90
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,11 @@ MODULE diag_data_mod
INTEGER, PARAMETER :: DIAG_SECONDS = 1, DIAG_MINUTES = 2, DIAG_HOURS = 3
INTEGER, PARAMETER :: DIAG_DAYS = 4, DIAG_MONTHS = 5, DIAG_YEARS = 6
INTEGER, PARAMETER :: MAX_SUBAXES = 10
INTEGER, PARAMETER :: NO_DOMAIN = 1 !< Use the FmsNetcdfFile_t fileobj
INTEGER, PARAMETER :: TWO_D_DOMAIN = 2 !< Use the FmsNetcdfDomainFile_t fileobj
INTEGER, PARAMETER :: UG_DOMAIN = 3 !< Use the FmsNetcdfUnstructuredDomainFile_t fileobj
INTEGER, PARAMETER :: DIRECTION_UP = 1 !< The axis points up if positive
INTEGER, PARAMETER :: DIRECTION_DOWN = -1 !< The axis points down if positive
INTEGER, PARAMETER :: GLO_REG_VAL = -999 !< Value used in the region specification of the diag_table
!! to indicate to use the full axis instead of a sub-axis
INTEGER, PARAMETER :: GLO_REG_VAL_ALT = -1 !< Alternate value used in the region specification of the
Expand Down
144 changes: 138 additions & 6 deletions diag_manager/fms_diag_axis_object.F90
Original file line number Diff line number Diff line change
Expand Up @@ -31,14 +31,17 @@ module fms_diag_axis_object_mod
use mpp_domains_mod, only: domain1d, domain2d, domainUG, mpp_get_compute_domain, CENTER, &
& mpp_get_compute_domain, NORTH, EAST
use platform_mod, only: r8_kind, r4_kind
use diag_data_mod, only: diag_atttype, max_axes
use diag_data_mod, only: diag_atttype, max_axes, NO_DOMAIN, TWO_D_DOMAIN, UG_DOMAIN, &
direction_down, direction_up
use mpp_mod, only: FATAL, mpp_error, uppercase
use fms2_io_mod, only: FmsNetcdfFile_t, FmsNetcdfDomainFile_t, FmsNetcdfUnstructuredDomainFile_t, &
& register_axis, register_field, register_variable_attribute, write_data
implicit none

PRIVATE

public :: diagAxis_t, set_subaxis, fms_diag_axis_init, fms_diag_axis_object_init, fms_diag_axis_object_end
public :: axis_obj
public :: diagAxis_t, set_subaxis, fms_diag_axis_init, fms_diag_axis_object_init, fms_diag_axis_object_end, &
& get_domain_and_domain_type, axis_obj
!> @}

!> @brief Type to hold the domain info for an axis
Expand Down Expand Up @@ -85,10 +88,13 @@ module fms_diag_axis_object_mod
CHARACTER(len=:), ALLOCATABLE, private :: long_name !< Long_name attribute of the axis
CHARACTER(len=1) , private :: cart_name !< Cartesian name "X", "Y", "Z", "T", "U", "N"
CLASS(*), ALLOCATABLE, private :: axis_data(:) !< Data of the axis
CHARACTER(len=:), ALLOCATABLE, private :: type_of_data !< The type of the axis_data ("float" or "double")
!< TO DO this can be a dlinked to avoid having limits
type(subaxis_t) , private :: subaxis(3) !< Array of subaxis
integer , private :: nsubaxis !< Number of subaxis
class(diagDomain_t),ALLOCATABLE, private :: axis_domain !< Domain
INTEGER , private :: type_of_domain !< The type of domain ("NO_DOMAIN", "TWO_D_DOMAIN",
!! or "UG_DOMAIN")
INTEGER , private :: length !< Global axis length
INTEGER , private :: direction !< Direction of the axis 0, 1, -1
INTEGER , private :: edges !< Axis ID for the previously defined "edges axis"
Expand All @@ -105,11 +111,10 @@ module fms_diag_axis_object_mod
PROCEDURE :: register => register_diag_axis_obj
PROCEDURE :: axis_length => get_axis_length
PROCEDURE :: set_subaxis
PROCEDURE :: write_axis_metadata
PROCEDURE :: write_axis_data

! TO DO:
! PROCEDURE :: write_axis_metadata
! PROCEDURE :: write_axis_data
! PROCEDURE :: get_fileobj_type_needed (use the domain to figure out what fms2 fileobj to use)
! Get/has/is subroutines as needed
END TYPE diagAxis_t

Expand Down Expand Up @@ -154,14 +159,17 @@ subroutine register_diag_axis_obj(obj, axis_name, axis_data, units, cart_name, l
type is (real(kind=r8_kind))
allocate(real(kind=r8_kind) :: obj%axis_data(size(axis_data)))
obj%axis_data = axis_data
obj%type_of_data = "double" !< This is what fms2_io expects in the register_field call
type is (real(kind=r4_kind))
allocate(real(kind=r4_kind) :: obj%axis_data(size(axis_data)))
obj%axis_data = axis_data
obj%type_of_data = "float" !< This is what fms2_io expects in the register_field call
class default
call mpp_error(FATAL, "The axis_data in your diag_axis_init call is not a supported type. &
& Currently only r4 and r8 data is supported.")
end select

obj%type_of_domain = NO_DOMAIN
if (present(Domain)) then
if (present(Domain2) .or. present(DomainU)) call mpp_error(FATAL, &
"The presence of Domain with any other domain type is prohibited. "//&
Expand All @@ -174,9 +182,11 @@ subroutine register_diag_axis_obj(obj, axis_name, axis_data, units, cart_name, l
"Check you diag_axis_init call for axis_name:"//trim(axis_name))
allocate(diagDomain2d_t :: obj%axis_domain)
call obj%axis_domain%set(Domain2=Domain2)
obj%type_of_domain = TWO_D_DOMAIN
else if (present(DomainU)) then
allocate(diagDomainUg_t :: obj%axis_domain)
call obj%axis_domain%set(DomainU=DomainU)
obj%type_of_domain = UG_DOMAIN
endif

obj%tile_count = 1
Expand All @@ -202,6 +212,96 @@ subroutine register_diag_axis_obj(obj, axis_name, axis_data, units, cart_name, l
obj%nsubaxis = 0
end subroutine register_diag_axis_obj

!> @brief Write the axis meta data to an open fileobj
subroutine write_axis_metadata(obj, fileobj, sub_axis_id)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this routine needs some more in-line documentation to make it easier for people to follow in the future. There are a lot of selects and cases to untangle, so general statements like "find the domain and write the correct axis type" might be useful

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This routine needs something to call it at the fmsDiagFileObject_type level. The axis should be added to the file type so this can be called using the fileobject for the file.

class(diagAxis_t), target, INTENT(IN) :: obj !< diag_axis obj
class(FmsNetcdfFile_t), INTENT(INOUT) :: fileobj !< Fms2_io fileobj to write the data to
integer, OPTIONAL, INTENT(IN) :: sub_axis_id !< ID of the sub_axis, if it exists

character(len=:), ALLOCATABLE :: axis_edges_name !< Name of the edges, if it exist
character(len=:), pointer :: axis_name !< Name of the axis
integer :: axis_length !< Size of the axis

if (present(sub_axis_id)) then
axis_name => obj%subaxis(sub_axis_id)%subaxis_name
axis_length = obj%subaxis(sub_axis_id)%ending_index - obj%subaxis(sub_axis_id)%starting_index + 1
else
axis_name => obj%axis_name
axis_length = obj%length
endif

!< Add the axis as a dimension in the netcdf file based on the type of axis_domain and the fileobj type
select type (fileobj)
type is (FmsNetcdfFile_t)
!< Here the axis is not domain decomposed (i.e z_axis)
call register_axis(fileobj, axis_name, axis_length)
type is (FmsNetcdfDomainFile_t)
select case (obj%type_of_domain)
case (NO_DOMAIN)
!< Here the fileobj is domain decomposed, but the axis is not
!! Domain decomposed fileobjs can have axis that are not domain decomposed (i.e "Z" axis)
call register_axis(fileobj, axis_name, axis_length)
case (TWO_D_DOMAIN)
!< Here the axis is domain decomposed
call register_axis(fileobj, axis_name, obj%cart_name, domain_position=obj%domain_position)
end select
type is (FmsNetcdfUnstructuredDomainFile_t)
select case (obj%type_of_domain)
case (NO_DOMAIN)
!< Here the fileobj is in the unstructured domain, but the axis is not
!< Unstructured domain fileobjs can have axis that are not domain decomposed (i.e "Z" axis)
call register_axis(fileobj, axis_name, axis_length)
case (UG_DOMAIN)
!< Here the axis is in a unstructured domain
call register_axis(fileobj, axis_name)
end select
end select

!< Add the axis as a variable and write its metada
call register_field(fileobj, axis_name, obj%type_of_data, (/axis_name/))
call register_variable_attribute(fileobj, axis_name, "longname", obj%long_name, &
str_len=len_trim(obj%long_name))

if (obj%cart_name .NE. "N") &
call register_variable_attribute(fileobj, axis_name, "axis", obj%cart_name, str_len=1)

if (trim(obj%units) .NE. "none") &
call register_variable_attribute(fileobj, axis_name, "units", obj%units, str_len=len_trim(obj%units))

select case (obj%direction)
case (direction_up)
call register_variable_attribute(fileobj, axis_name, "positive", "up", str_len=2)
case (direction_down)
call register_variable_attribute(fileobj, axis_name, "positive", "down", str_len=4)
end select

if (obj%edges > 0) then
axis_edges_name = axis_obj(obj%edges)%axis_name
call register_variable_attribute(fileobj, axis_name, "edges", axis_edges_name, &
str_len=len_trim(axis_edges_name))
endif

end subroutine write_axis_metadata

!> @brief Write the axis data to an open fileobj
subroutine write_axis_data(obj, fileobj, sub_axis_id)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This also needs a routine to call it from the file object level.

class(diagAxis_t), INTENT(IN) :: obj !< diag_axis obj
class(FmsNetcdfFile_t), INTENT(INOUT) :: fileobj !< Fms2_io fileobj to write the data to
integer, OPTIONAL, INTENT(IN) :: sub_axis_id !< ID of the sub_axis, if it exists

integer :: i !< Starting index of a sub_axis
integer :: j !< Ending index of a sub_axis

if (present(sub_axis_id)) then
i = obj%subaxis(sub_axis_id)%starting_index
j = obj%subaxis(sub_axis_id)%ending_index

call write_data(fileobj, obj%subaxis(sub_axis_id)%subaxis_name, obj%axis_data(i:j))
else
call write_data(fileobj, obj%axis_name, obj%axis_data)
endif
end subroutine write_axis_data

!> @brief Get the length of the axis
!> @return axis length
function get_axis_length(obj) &
Expand Down Expand Up @@ -384,6 +484,38 @@ subroutine check_if_valid_edges(edges)
call mpp_error(FATAL, "diag_axit_init: The edge axis has not been defined. "&
"Call diag_axis_init for the edge axis first")
end subroutine check_if_valid_edges

!> @brief Loop through a variable's axis_id to determine and return the domain type and domain to use
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still think this description is lacking. The \brief describes what this subroutine does, but there is no explanation for it's use. I think with the next update, you should have a \description instead of \brief that explains why you use this. In the larger plan, this routine is necessary. If you're just looking at this module, you will scratch your head and wonder why this exists.

subroutine get_domain_and_domain_type(axis_id, domain_type, domain, var_name)
integer, INTENT(IN) :: axis_id(:) !< Array of axis ids
integer, INTENT(OUT) :: domain_type !< fileobj_type to use
CLASS(diagDomain_t), POINTER, INTENT(OUT) :: domain !< Domain
character(len=*), INTENT(IN) :: var_name !< Name of the variable (for error messages)

integer :: i !< For do loops
integer :: j !< axis_id(i) (for less typing)

domain_type = NO_DOMAIN
domain => null()

do i = 1, size(axis_id)
j = axis_id(i)
!< Check that all the axis are in the same domain
if (domain_type .ne. axis_obj(j)%type_of_domain) then
!< If they are different domains, one of them can be NO_DOMAIN
!! i.e a variable can have axis that are domain decomposed (x,y) and an axis that isn't (z)
if (domain_type .eq. NO_DOMAIN .or. axis_obj(j)%type_of_domain .eq. NO_DOMAIN ) then
!< Update the domain_type and domain, if needed
if (axis_obj(j)%type_of_domain .eq. TWO_D_DOMAIN .or. axis_obj(j)%type_of_domain .eq. UG_DOMAIN) then
domain_type = axis_obj(j)%type_of_domain
domain => axis_obj(j)%axis_domain
endif
else
call mpp_error(FATAL, "The variable:"//trim(var_name)//" has axis that are not in the same domain")
endif
endif
enddo
end subroutine get_domain_and_domain_type
end module fms_diag_axis_object_mod
!> @}
! close documentation grouping
11 changes: 11 additions & 0 deletions diag_manager/fms_diag_object.F90
Original file line number Diff line number Diff line change
Expand Up @@ -227,6 +227,17 @@ subroutine fms_register_diag_field_obj &
! return
! endif

!> TO DO: Add all the info from the diag_axis obj
!! axes will need to be changed to optional, so this subroutine can be used for both scalar and array fields
!! the domain_type and domain will be need to added to the dobj
! if (present(axes))
! dobj%axes => axes ! or something
! call get_domain_and_domain_type(dobj%axes, dobj%domain_type, dobj%domain, dobj%varname)
!! Send all the axes_info to the diag_files
! else
! dobj%domain_type = NO_DOMAIN
! endif

!> get the optional arguments if included and the diagnostic is in the diag table
if (present(longname)) then
allocate(character(len=len(longname)) :: dobj%longname)
Expand Down