Merge pull request #3212 from GEOS-ESM/feature/fixes-#3175

Feature/fixes #3175
GEOS-ESM · Dec 4, 2024 · 6f84fdb · 6f84fdb
2 parents 58ddb04 + 9a2fee5
commit 6f84fdb
Show file tree

Hide file tree

Showing 3 changed files with 135 additions and 23 deletions.
diff --git 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 &#177;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
@@ -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
@@ -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)