diff --git a/CHANGELOG.md b/CHANGELOG.md index 56af8fffc55d..b8234715b962 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Allow update offsets of ±timestep in ExtData2G - Minor revision (and generalization) of grid-def for GSI purposes - Trajectory sampler: fix a bug when group_name does not exist in netCDF file and a bug that omitted the first time point +- Allow lat-lon grid factory to detect and use CF compliant lat-lon bounds in a file when making a grid - PFIO/Variable class, new procedures to retrieve string/reals/int attributes from a variable ### Changed @@ -45,6 +46,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed - Fixed issue of some Baselibs builds appearing to support zstandard. This is not possible due to Baselibs building HDF5 and netCDF as static libraries +- Fixed a bug where the periodicity around the earth of the lat-lon grid was not being set properly when grid did not span from pole to pole ### Removed diff --git a/base/MAPL_LatLonGridFactory.F90 b/base/MAPL_LatLonGridFactory.F90 index fbbbfe3a41e7..b7d076e5a2a3 100644 --- a/base/MAPL_LatLonGridFactory.F90 +++ b/base/MAPL_LatLonGridFactory.F90 @@ -56,6 +56,8 @@ module MAPL_LatLonGridFactoryMod integer :: px, py logical :: is_halo_initialized = .false. logical :: periodic = .true. + character(len=:), allocatable :: lon_bounds_name + character(len=:), allocatable :: lat_bounds_name contains procedure :: make_new_grid procedure :: create_basic_grid @@ -218,10 +220,16 @@ function create_basic_grid(this, unusable, rc) result(grid) integer, optional, intent(out) :: rc integer :: status + type(ESMF_PoleKind_Flag) :: polekindflag(2) _UNUSED_DUMMY(unusable) if (this%periodic) then + if (this%pole == "XY") then + polekindflag = ESMF_POLEKIND_NONE + else + polekindflag = ESMF_POLEKIND_MONOPOLE + end if grid = ESMF_GridCreate1PeriDim( & & name = this%grid_name, & & countsPerDEDim1=this%ims, & @@ -232,6 +240,7 @@ function create_basic_grid(this, unusable, rc) result(grid) & coordDep1=[1,2], & & coordDep2=[1,2], & & coordSys=ESMF_COORDSYS_SPH_RAD, & + & polekindflag=polekindflag, & & rc=status) _VERIFY(status) else @@ -673,7 +682,7 @@ subroutine initialize_from_file_metadata(this, file_metadata, unusable, force_fi integer :: i_min, i_max real(kind=REAL64) :: d_lat, d_lat_temp, extrap_lat - logical :: is_valid, use_file_coords, compute_lons, compute_lats + logical :: is_valid, use_file_coords, compute_lons, compute_lats, has_bnds _UNUSED_DUMMY(unusable) @@ -759,6 +768,11 @@ subroutine initialize_from_file_metadata(this, file_metadata, unusable, force_fi where(this%lon_centers > 180) this%lon_centers=this%lon_centers-360 end if + has_bnds = coordinate_has_bounds(file_metadata, lon_name, _RC) + if (has_bnds) then + this%lon_bounds_name = get_coordinate_bounds_name(file_metadata, lon_name, _RC) + this%lon_corners = get_coordinate_bounds(file_metadata, lon_name, _RC) + end if v => file_metadata%get_coordinate_variable(lat_name, rc=status) _VERIFY(status) @@ -773,6 +787,12 @@ subroutine initialize_from_file_metadata(this, file_metadata, unusable, force_fi _FAIL('unsupported type of data; must be REAL32 or REAL64') end select + has_bnds = coordinate_has_bounds(file_metadata, lat_name, _RC) + if (has_bnds) then + this%lat_bounds_name = get_coordinate_bounds_name(file_metadata, lat_name, _RC) + this%lat_corners = get_coordinate_bounds(file_metadata, lat_name, _RC) + end if + ! Check: is this a "mis-specified" pole-centered grid? if (size(this%lat_centers) >= 4) then @@ -804,14 +824,14 @@ subroutine initialize_from_file_metadata(this, file_metadata, unusable, force_fi end if end if - ! Corners are the midpoints of centers (and extrapolated at the ! poles for lats.) - allocate(this%lon_corners(im+1), this%lat_corners(jm+1)) - - this%lon_corners(1) = (this%lon_centers(im) + this%lon_centers(1))/2 - 180 - this%lon_corners(2:im) = (this%lon_centers(1:im-1) + this%lon_centers(2:im))/2 - this%lon_corners(im+1) = (this%lon_centers(im) + this%lon_centers(1))/2 + 180 + if (.not. allocated(this%lon_corners)) then + allocate(this%lon_corners(im+1)) + this%lon_corners(1) = (this%lon_centers(im) + this%lon_centers(1))/2 - 180 + this%lon_corners(2:im) = (this%lon_centers(1:im-1) + this%lon_centers(2:im))/2 + this%lon_corners(im+1) = (this%lon_centers(im) + this%lon_centers(1))/2 + 180 + end if ! This section about pole/dateline is probably not needed in file data case. if (abs(this%lon_centers(1) + 180) < 1000*epsilon(1.0)) then @@ -826,10 +846,13 @@ subroutine initialize_from_file_metadata(this, file_metadata, unusable, force_fi this%dateline = 'XY' this%lon_range = RealMinMax(this%lon_centers(1), this%lon_centers(jm)) end if - - this%lat_corners(1) = this%lat_centers(1) - (this%lat_centers(2)-this%lat_centers(1))/2 - this%lat_corners(2:jm) = (this%lat_centers(1:jm-1) + this%lat_centers(2:jm))/2 - this%lat_corners(jm+1) = this%lat_centers(jm) - (this%lat_centers(jm-1)-this%lat_centers(jm))/2 + + if (.not. allocated(this%lat_corners)) then + allocate(this%lat_corners(jm+1)) + this%lat_corners(1) = this%lat_centers(1) - (this%lat_centers(2)-this%lat_centers(1))/2 + this%lat_corners(2:jm) = (this%lat_centers(1:jm-1) + this%lat_centers(2:jm))/2 + this%lat_corners(jm+1) = this%lat_centers(jm) - (this%lat_centers(jm-1)-this%lat_centers(jm))/2 + end if if (abs(this%lat_centers(1) + 90) < 1000*epsilon(1.0)) then this%pole = 'PC' @@ -1139,7 +1162,6 @@ subroutine check_and_fill_consistency(this, unusable, rc) ! Check regional vs global if (this%pole == 'XY') then ! regional - this%periodic = .false. _ASSERT(this%lat_range%min /= MAPL_UNDEFINED_REAL, 'uninitialized min for lat_range') _ASSERT(this%lat_range%max /= MAPL_UNDEFINED_REAL, 'uninitialized min for lat_range') else ! global @@ -1849,9 +1871,16 @@ function get_file_format_vars(this) result(vars) class (LatLonGridFactory), intent(inout) :: this character(len=:), allocatable :: vars + integer :: i _UNUSED_DUMMY(this) vars = 'lon,lat' + if (allocated(this%lon_bounds_name)) then + vars = vars // ',' // this%lon_bounds_name + end if + if (allocated(this%lat_bounds_name)) then + vars = vars // ',' // this%lat_bounds_name + end if end function get_file_format_vars @@ -1928,5 +1957,85 @@ function generate_file_reference3D(this,fpointer,metaData) result(ref) _UNUSED_DUMMY(metaData) end function generate_file_reference3D + function coordinate_has_bounds(metadata, coord_name, rc) result(has_bounds) + logical :: has_bounds + type(FileMetadata), intent(in) :: metadata + character(len=*), intent(in) :: coord_name + integer, optional, intent(out) :: rc + + type(Variable), pointer :: var + integer :: status + + var => metadata%get_variable(coord_name, _RC) + has_bounds = var%is_attribute_present("bounds") + + _RETURN(_SUCCESS) + end function + + function get_coordinate_bounds_name(metadata, coord_name, rc) result(coord_bounds_name) + character(len=:), allocatable :: coord_bounds_name + type(FileMetadata), intent(in) :: metadata + character(len=*), intent(in) :: coord_name + integer, optional, intent(out) :: rc + + type(Variable), pointer :: var + type(Attribute), pointer :: attr + integer :: status + class(*), pointer :: attr_val + + var => metadata%get_variable(coord_name, _RC) + attr => var%get_attribute("bounds", _RC) + attr_val => attr%get_value() + select type(attr_val) + type is(character(*)) + coord_bounds_name = attr_val + class default + _FAIL('coordinate bounds must be a string') + end select + _RETURN(_SUCCESS) + end function + + function get_coordinate_bounds(metadata, coord_name, rc) result(coord_bounds) + real(kind=REAL64), allocatable :: coord_bounds(:) + type(FileMetadata), intent(in) :: metadata + character(len=*), intent(in) :: coord_name + integer, optional, intent(out) :: rc + + type(Variable), pointer :: var + type(Attribute), pointer :: attr + integer :: status, im, i + class(*), pointer :: attr_val + character(len=:), allocatable :: bnds_name, source_file + real(kind=REAL64), allocatable :: file_bounds(:,:) + type(NetCDF4_FileFormatter) :: file_formatter + + + var => metadata%get_variable(coord_name, _RC) + attr => var%get_attribute("bounds", _RC) + attr_val => attr%get_value() + select type(attr_val) + type is(character(*)) + bnds_name = attr_val + class default + _FAIL('coordinate bounds must be a string') + end select + im = metadata%get_dimension(coord_name, _RC) + allocate(coord_bounds(im+1), _STAT) + allocate(file_bounds(2,im), _STAT) + source_file = metadata%get_source_file() + + call file_formatter%open(source_file, PFIO_READ, _RC) + call file_formatter%get_var(bnds_name, file_bounds, _RC) + call file_formatter%close(_RC) + do i=1,im-1 + _ASSERT(file_bounds(2,i)==file_bounds(1,i+1), "Bounds are not contiguous in file") + enddo + do i=1,im + coord_bounds(i) = file_bounds(1,i) + coord_bounds(i+1) = file_bounds(2,i) + enddo + + _RETURN(_SUCCESS) + end function end module MAPL_LatLonGridFactoryMod diff --git a/griddedio/FieldBundleRead.F90 b/griddedio/FieldBundleRead.F90 index 6f0bd2b09c65..352dd414c330 100644 --- a/griddedio/FieldBundleRead.F90 +++ b/griddedio/FieldBundleRead.F90 @@ -48,8 +48,8 @@ subroutine MAPL_create_bundle_from_metdata_id(bundle,metadata_id,file_name,only_ type(StringVariableMap), pointer :: variables type(Variable), pointer :: this_variable type(StringVariableMapIterator) :: var_iter - character(len=:), pointer :: var_name,dim_name - character(len=:), allocatable :: lev_name + character(len=:), pointer :: var_name_ptr,dim_name + character(len=:), allocatable :: lev_name,var_name type(ESMF_Field) :: field type (StringVector), pointer :: dimensions type (StringVectorIterator) :: dim_iter @@ -71,14 +71,15 @@ subroutine MAPL_create_bundle_from_metdata_id(bundle,metadata_id,file_name,only_ factory => get_factory(file_grid,rc=status) _VERIFY(status) grid_vars = factory%get_file_format_vars() - exclude_vars = grid_vars//",lev,time,lons,lats" + exclude_vars = ","//grid_vars//",lev,time,time_bnds," if (has_vertical_level) lev_size = metadata%get_dimension(trim(lev_name)) variables => metadata%get_variables() var_iter = variables%begin() do while (var_iter /= variables%end()) var_has_levels = .false. - var_name => var_iter%key() + var_name_ptr => var_iter%key() + var_name = ","//var_name_ptr//"," this_variable => var_iter%value() if (has_vertical_level) then @@ -91,20 +92,20 @@ subroutine MAPL_create_bundle_from_metdata_id(bundle,metadata_id,file_name,only_ enddo end if - if (index(','//trim(exclude_vars)//',',','//trim(var_name)//',') > 0) then + if (index(trim(exclude_vars),trim(var_name)) > 0) then call var_iter%next() cycle end if create_variable = .true. if (present(only_vars)) then - if (index(','//trim(only_vars)//',',','//trim(var_name)//',') < 1) create_variable = .false. + if (index(','//trim(only_vars)//',',trim(var_name)) < 1) create_variable = .false. end if if (create_variable) then if(var_has_levels) then if (grid_size(3) == lev_size) then location=MAPL_VLocationCenter dims = MAPL_DimsHorzVert - field= ESMF_FieldCreate(grid,name=trim(var_name),typekind=ESMF_TYPEKIND_R4, & + field= ESMF_FieldCreate(grid,name=trim(var_name_ptr),typekind=ESMF_TYPEKIND_R4, & ungriddedUbound=[grid_size(3)],ungriddedLBound=[1], rc=status) block real, pointer :: ptr3d(:,:,:) @@ -114,7 +115,7 @@ subroutine MAPL_create_bundle_from_metdata_id(bundle,metadata_id,file_name,only_ else if (grid_size(3)+1 == lev_size) then location=MAPL_VLocationEdge dims = MAPL_DimsHorzVert - field= ESMF_FieldCreate(grid,name=trim(var_name),typekind=ESMF_TYPEKIND_R4, & + field= ESMF_FieldCreate(grid,name=trim(var_name_ptr),typekind=ESMF_TYPEKIND_R4, & ungriddedUbound=[grid_size(3)],ungriddedLBound=[0], rc=status) block real, pointer :: ptr3d(:,:,:) @@ -125,7 +126,7 @@ subroutine MAPL_create_bundle_from_metdata_id(bundle,metadata_id,file_name,only_ else location=MAPL_VLocationNone dims = MAPL_DimsHorzOnly - field= ESMF_FieldCreate(grid,name=trim(var_name),typekind=ESMF_TYPEKIND_R4, & + field= ESMF_FieldCreate(grid,name=trim(var_name_ptr),typekind=ESMF_TYPEKIND_R4, & rc=status) block real, pointer :: ptr2d(:,:) @@ -137,8 +138,8 @@ subroutine MAPL_create_bundle_from_metdata_id(bundle,metadata_id,file_name,only_ _VERIFY(status) call ESMF_AttributeSet(field,name='VLOCATION',value=location,rc=status) _VERIFY(status) - units = metadata%get_var_attr_string(var_name,'units',_RC) - long_name = metadata%get_var_attr_string(var_name,'long_name',_RC) + units = metadata%get_var_attr_string(var_name_ptr,'units',_RC) + long_name = metadata%get_var_attr_string(var_name_ptr,'long_name',_RC) call ESMF_AttributeSet(field,name='UNITS',value=units,rc=status) _VERIFY(status) call ESMF_AttributeSet(field,name='LONG_NAME',value=long_name,rc=status)