Feature/wjiang/cs to loc

CHANGELOG.md

@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Added test case for regridding from CubedSphere to Loctream
 - Added regression test for Regrid\_Util.x
 - External pfio server GridComp and ctest: new `mapl_PfioServerGridComp_mod` provides
   an ESMF GridComp whose `run` phase creates and starts an `MpiServer` or

gridcomps/componentDriverGridComp/componentDriverGridComp.F90

@@ -200,6 +200,36 @@ subroutine initialize_internal_state(internal_state, support, hconfig, rc)
       type(ESMF_HConfig), intent(in) :: hconfig
       integer, optional, intent(out) :: rc
 
+      integer :: status
+      type(ESMF_Field) :: field
+      type(ESMF_GeomType_Flag) :: geomtype
+
+      ! The internal-state fields (rand, grid_lons, grid_lats, quarter_grid)
+      ! are built on this component's own geometry. That geometry is normally
+      ! an ESMF_Grid (rank-2 horizontal fields), but it may also be a
+      ! non-gridded geometry such as ESMF_LocStream (rank-1 fields). Dispatch
+      ! on the actual geom type of one of these fields to pick the correct
+      ! rank-specific implementation.
+      call ESMF_StateGet(internal_state, 'rand', field, _RC)
+      call ESMF_FieldGet(field, geomtype=geomtype, _RC)
+
+      if (geomtype == ESMF_GEOMTYPE_GRID) then
+         call initialize_internal_state_grid(internal_state, hconfig, _RC)
+      else if (geomtype == ESMF_GEOMTYPE_LOCSTREAM) then
+         call initialize_internal_state_locstream(internal_state, hconfig, _RC)
+      else
+         _FAIL('componentDriverGridComp: unsupported geometry type for internal state')
+      end if
+
+      _RETURN(_SUCCESS)
+
+   end subroutine initialize_internal_state
+
+   subroutine initialize_internal_state_grid(internal_state, hconfig, rc)
+      type(ESMF_State), intent(inout) :: internal_state
+      type(ESMF_HConfig), intent(in) :: hconfig
+      integer, optional, intent(out) :: rc
+
       real, pointer :: ptr_2d(:, :)
       real(kind=ESMF_KIND_R8), pointer :: coords(:, :)
       integer :: status, seed_size, mypet, i, j
@@ -253,7 +283,61 @@ subroutine initialize_internal_state(internal_state, support, hconfig, rc)
 
       _RETURN(_SUCCESS)
 
-   end subroutine initialize_internal_state
+   end subroutine initialize_internal_state_grid
+
+   subroutine initialize_internal_state_locstream(internal_state, hconfig, rc)
+      type(ESMF_State), intent(inout) :: internal_state
+      type(ESMF_HConfig), intent(in) :: hconfig
+      integer, optional, intent(out) :: rc
+
+      real, pointer :: ptr_1d(:)
+      real(kind=ESMF_KIND_R8), pointer :: coords(:)
+      integer :: status, seed_size, mypet, i
+      integer, allocatable :: seeds(:)
+      type(ESMF_Field) :: field
+      type(ESMF_LocStream) :: locstream
+      type(ESMF_VM) :: vm
+      logical :: is_present
+      real :: quarter_grid_fac1, quarter_grid_fac2
+
+      ! rand
+      call MAPL_StateGetPointer(internal_state, ptr_1d, 'rand', _RC)
+      call random_seed(size=seed_size)
+      allocate(seeds(seed_size))
+      call ESMF_VMGetCurrent(vm, _RC)
+      call ESMF_VMGet(vm, localPet=mypet, _RC)
+      seeds = mypet
+      call random_seed(put=seeds)
+      call random_number(ptr_1d)
+      ! lons and lats
+      call MAPL_StateGetPointer(internal_state, ptr_1d, 'grid_lons', _RC)
+      call ESMF_StateGet(internal_state, 'grid_lons', field, _RC)
+      call ESMF_FieldGet(field, locstream=locstream, _RC)
+      call ESMF_LocStreamGetKey(locstream, keyName='ESMF:Lon', farray=coords, _RC)
+      ptr_1d = real(coords, kind=kind(ptr_1d))
+      call MAPL_StateGetPointer(internal_state, ptr_1d, 'grid_lats', _RC)
+      call ESMF_LocStreamGetKey(locstream, keyName='ESMF:Lat', farray=coords, _RC)
+      ptr_1d = real(coords, kind=kind(ptr_1d))
+
+      quarter_grid_fac1 = 1.0
+      quarter_grid_fac2 = 2.0
+      is_present = ESMF_HConfigIsDefined(hconfig, keyString='quarter_grid_fac1', _RC)
+      if (is_present) then
+         quarter_grid_fac1 = ESMF_HConfigAsR4(hconfig, keyString='quarter_grid_fac1', _RC)
+      end if
+      is_present = ESMF_HConfigIsDefined(hconfig, keyString='quarter_grid_fac2', _RC)
+      if (is_present) then
+         quarter_grid_fac2 = ESMF_HConfigAsR4(hconfig, keyString='quarter_grid_fac2', _RC)
+      end if
+      call MAPL_StateGetPointer(internal_state, ptr_1d, 'quarter_grid', _RC)
+      ptr_1d = quarter_grid_fac2
+      do i = 1, size(ptr_1d), 2
+         ptr_1d(i) = quarter_grid_fac1
+      end do
+
+      _RETURN(_SUCCESS)
+
+   end subroutine initialize_internal_state_locstream
 
    subroutine update_internal_state(internal_state, current_time, support, rc)
       type(ESMF_State), intent(inout) :: internal_state
@@ -263,9 +347,28 @@ subroutine update_internal_state(internal_state, current_time, support, rc)
 
       integer :: status
       real, pointer :: ptr_2d(:, :)
-
-      call MAPL_StateGetPointer(internal_state, ptr_2d, 'time_interval', _RC)
-      ptr_2d = support%tFunc%evaluate_time(current_time, _RC)
+      real, pointer :: ptr_1d(:)
+      type(ESMF_Field) :: field
+      type(ESMF_GeomType_Flag) :: geomtype
+      real :: time_value
+
+      ! The 'time_interval' internal field is built on this component's own
+      ! geometry, which may be an ESMF_Grid (rank-2) or a non-gridded
+      ! geometry such as ESMF_LocStream (rank-1).
+      call ESMF_StateGet(internal_state, 'time_interval', field, _RC)
+      call ESMF_FieldGet(field, geomtype=geomtype, _RC)
+
+      time_value = support%tFunc%evaluate_time(current_time, _RC)
+
+      if (geomtype == ESMF_GEOMTYPE_GRID) then
+         call MAPL_StateGetPointer(internal_state, ptr_2d, 'time_interval', _RC)
+         ptr_2d = time_value
+      else if (geomtype == ESMF_GEOMTYPE_LOCSTREAM) then
+         call MAPL_StateGetPointer(internal_state, ptr_1d, 'time_interval', _RC)
+         ptr_1d = time_value
+      else
+         _FAIL('componentDriverGridComp: unsupported geometry type for internal state')
+      end if
 
       _RETURN(_SUCCESS)
 

infrastructure/esmf/FieldPointerUtilities.F90

@@ -428,30 +428,41 @@ subroutine clone(x, y, name, rc)
 
       character(len=*), parameter :: CLONE_TAG = '_clone'
       !type(ESMF_ArraySpec) :: arrayspec
+      type(ESMF_Geom) :: geom
       type(ESMF_Grid) :: grid
+      type(ESMF_GeomType_Flag) :: geomtype
       type(ESMF_StaggerLoc) :: staggerloc
       integer, allocatable :: gridToFieldMap(:)
       integer, allocatable :: ungriddedLBound(:)
       integer, allocatable :: ungriddedUBound(:)
       type(ESMF_TypeKind_Flag) :: tk
       character(len=ESMF_MAXSTR) :: clone_name
       integer :: status
-      integer :: field_rank, grid_rank,ungrid_size
+      integer :: field_rank, geom_rank, ungrid_size
+      logical :: is_grid
       type(ESMF_Index_Flag) :: index_flag
       real(kind=ESMF_KIND_R4), pointer      :: VR4_1D(:), VR4_2D(:,:), VR4_3D(:,:,:), VR4_4D(:,:,:,:)
       real(kind=ESMF_KIND_R8), pointer      :: VR8_1D(:), VR8_2D(:,:), VR8_3D(:,:,:), VR8_4D(:,:,:,:)
       integer, allocatable :: lc(:)
       type(ESMF_Info) :: x_info, y_info
 
-      call ESMF_FieldGet(x,grid=grid,rank=field_rank,_RC)
+      call ESMF_FieldGet(x, geom=geom, rank=field_rank, _RC)
       lc = get_local_element_count(x,_RC)
-      call ESMF_GridGet(grid,dimCount=grid_rank,indexFlag=index_flag,_RC)
-      ungrid_size = field_rank-grid_rank
-      allocate(gridToFieldMap(grid_rank))
+      call ESMF_GeomGet(geom, geomtype=geomtype, dimCount=geom_rank, indexFlag=index_flag, _RC)
+      is_grid = (geomtype == ESMF_GEOMTYPE_GRID)
+      if (is_grid) call ESMF_GeomGet(geom, grid=grid, _RC)
+      ungrid_size = field_rank-geom_rank
+      allocate(gridToFieldMap(geom_rank))
       allocate(ungriddedLBound(ungrid_size),ungriddedUBound(ungrid_size))
-      call ESMF_FieldGet(x, typekind=tk, name=clone_name, &
-         staggerloc=staggerloc, gridToFieldMap=gridToFieldMap, &
-         ungriddedLBound=ungriddedLBound, ungriddedUBound=ungriddedUBound,  _RC)
+      if (is_grid) then
+         call ESMF_FieldGet(x, typekind=tk, name=clone_name, &
+            staggerloc=staggerloc, gridToFieldMap=gridToFieldMap, &
+            ungriddedLBound=ungriddedLBound, ungriddedUBound=ungriddedUBound,  _RC)
+      else
+         call ESMF_FieldGet(x, typekind=tk, name=clone_name, &
+            gridToFieldMap=gridToFieldMap, &
+            ungriddedLBound=ungriddedLBound, ungriddedUBound=ungriddedUBound,  _RC)
+      end if
 
       if (present(name)) then
          clone_name = name
@@ -460,37 +471,73 @@ subroutine clone(x, y, name, rc)
       end if
 
       if (index_flag == ESMF_INDEX_USER) then
-         if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 1) then
-            allocate(VR4_1d(lc(1)),_STAT)
-            y = ESMF_FieldCreate(grid,VR4_1d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
-         else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 1) then
-            allocate(VR8_1d(lc(1)),_STAT)
-            y = ESMF_FieldCreate(grid,VR8_1d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
-         else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 2) then
-            allocate(VR4_2d(lc(1),lc(2)),_STAT)
-            y = ESMF_FieldCreate(grid,VR4_2d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
-         else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 2) then
-            allocate(VR8_2d(lc(1),lc(2)),_STAT)
-            y = ESMF_FieldCreate(grid,VR8_2d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
-         else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 3) then
-            allocate(VR4_3d(lc(1),lc(2),lc(3)),_STAT)
-            y = ESMF_FieldCreate(grid,VR4_3d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
-         else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 3) then
-            allocate(VR8_3d(lc(1),lc(2),lc(3)),_STAT)
-            y = ESMF_FieldCreate(grid,VR8_3d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
-         else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 4) then
-            allocate(VR4_4d(lc(1),lc(2),lc(3),lc(4)),_STAT)
-            y = ESMF_FieldCreate(grid,VR4_4d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
-         else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 4) then
-            allocate(VR8_4d(lc(1),lc(2),lc(3),lc(4)),_STAT)
-            y = ESMF_FieldCreate(grid,VR8_4d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+         if (is_grid) then
+            if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 1) then
+               allocate(VR4_1d(lc(1)),_STAT)
+               y = ESMF_FieldCreate(grid,VR4_1d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 1) then
+               allocate(VR8_1d(lc(1)),_STAT)
+               y = ESMF_FieldCreate(grid,VR8_1d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 2) then
+               allocate(VR4_2d(lc(1),lc(2)),_STAT)
+               y = ESMF_FieldCreate(grid,VR4_2d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 2) then
+               allocate(VR8_2d(lc(1),lc(2)),_STAT)
+               y = ESMF_FieldCreate(grid,VR8_2d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 3) then
+               allocate(VR4_3d(lc(1),lc(2),lc(3)),_STAT)
+               y = ESMF_FieldCreate(grid,VR4_3d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 3) then
+               allocate(VR8_3d(lc(1),lc(2),lc(3)),_STAT)
+               y = ESMF_FieldCreate(grid,VR8_3d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 4) then
+               allocate(VR4_4d(lc(1),lc(2),lc(3),lc(4)),_STAT)
+               y = ESMF_FieldCreate(grid,VR4_4d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 4) then
+               allocate(VR8_4d(lc(1),lc(2),lc(3),lc(4)),_STAT)
+               y = ESMF_FieldCreate(grid,VR8_4d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else
+               _FAIL( 'unsupported typekind+field_rank')
+            end if
          else
-            _FAIL( 'unsupported typekind+field_rank')
+            if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 1) then
+               allocate(VR4_1d(lc(1)),_STAT)
+               y = ESMF_FieldCreate(geom,VR4_1d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 1) then
+               allocate(VR8_1d(lc(1)),_STAT)
+               y = ESMF_FieldCreate(geom,VR8_1d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 2) then
+               allocate(VR4_2d(lc(1),lc(2)),_STAT)
+               y = ESMF_FieldCreate(geom,VR4_2d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 2) then
+               allocate(VR8_2d(lc(1),lc(2)),_STAT)
+               y = ESMF_FieldCreate(geom,VR8_2d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 3) then
+               allocate(VR4_3d(lc(1),lc(2),lc(3)),_STAT)
+               y = ESMF_FieldCreate(geom,VR4_3d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 3) then
+               allocate(VR8_3d(lc(1),lc(2),lc(3)),_STAT)
+               y = ESMF_FieldCreate(geom,VR8_3d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R4 .and. field_rank == 4) then
+               allocate(VR4_4d(lc(1),lc(2),lc(3),lc(4)),_STAT)
+               y = ESMF_FieldCreate(geom,VR4_4d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else if (tk == ESMF_TYPEKIND_R8 .and. field_rank == 4) then
+               allocate(VR8_4d(lc(1),lc(2),lc(3),lc(4)),_STAT)
+               y = ESMF_FieldCreate(geom,VR8_4d,gridToFieldMap=gridToFieldMap,name=clone_name,_RC)
+            else
+               _FAIL( 'unsupported typekind+field_rank')
+            end if
          end if
       else
-         y = ESMF_FieldCreate(grid, tk, staggerloc=staggerloc, &
-            gridToFieldMap=gridToFieldMap, ungriddedLBound=ungriddedLBound, &
-            ungriddedUBound=ungriddedUBound, name=clone_name, _RC)
+         if (is_grid) then
+            y = ESMF_FieldCreate(grid, tk, staggerloc=staggerloc, &
+               gridToFieldMap=gridToFieldMap, ungriddedLBound=ungriddedLBound, &
+               ungriddedUBound=ungriddedUBound, name=clone_name, _RC)
+         else
+            y = ESMF_FieldCreate(geom, tk, &
+               gridToFieldMap=gridToFieldMap, ungriddedLBound=ungriddedLBound, &
+               ungriddedUBound=ungriddedUBound, name=clone_name, _RC)
+         end if
       end if
 
       ! clone metadata

infrastructure/field_bundle/FieldBundleGet.F90

@@ -141,13 +141,19 @@ subroutine get_geom(fieldBundle, geom, rc)
          integer :: status
          type(ESMF_GeomType_Flag) :: geomtype
          type(ESMF_Grid) :: grid
+         type(ESMF_LocStream) :: locstream
 
          call ESMF_FieldBundleGet(fieldBundle, geomtype=geomtype, _RC)
          if (geomtype == ESMF_GEOMTYPE_GRID) then
             call ESMF_FieldBundleGet(fieldBundle, grid=grid, _RC)
             ! probable memory leak
             geom = ESMF_GeomCreate(grid=grid, _RC)
             _RETURN(_SUCCESS)
+         else if (geomtype == ESMF_GEOMTYPE_LOCSTREAM) then
+            call ESMF_FieldBundleGet(fieldBundle, locstream=locstream, _RC)
+            ! probable memory leak
+            geom = ESMF_GeomCreate(locstream=locstream, _RC)
+            _RETURN(_SUCCESS)
          end if
 
          _FAIL('unsupported geomtype; needs simple extension')

infrastructure/field_bundle/FieldBundleSet.F90

@@ -80,6 +80,7 @@ subroutine bundle_set(fieldBundle, unusable, &
       type(ESMF_GeomType_Flag) :: geomtype
       type(ESMF_Info) :: bundle_info
       type(ESMF_Grid) :: grid
+      type(ESMF_LocStream) :: locstream
       integer :: i
       type(ESMF_Field), allocatable :: fieldList(:)
       logical, allocatable :: has_geom
@@ -96,6 +97,15 @@ subroutine bundle_set(fieldBundle, unusable, &
             call FieldBundleReset(fieldBundle)
             call ESMF_FieldBundleSet(fieldBundle, grid=grid, _RC)
 
+            call FieldBundleGet(fieldBundle, fieldList=fieldList, _RC)
+            do i = 1, size(fieldList)
+               call MAPL_FieldSet(fieldList(i), geom=geom, _RC)
+            end do
+         else if (geomtype == ESMF_GEOMTYPE_LOCSTREAM) then
+            call ESMF_GeomGet(geom, locstream=locstream, _RC)
+            call FieldBundleReset(fieldBundle)
+            call ESMF_FieldBundleSet(fieldBundle, locstream=locstream, _RC)
+
             call FieldBundleGet(fieldBundle, fieldList=fieldList, _RC)
             do i = 1, size(fieldList)
                call MAPL_FieldSet(fieldList(i), geom=geom, _RC)

infrastructure/geom/GeomManager/initialize.F90

@@ -9,17 +9,20 @@
    module subroutine initialize(this)
       use mapl_LatLonGeomFactory_mod
       use mapl_CubedSphereGeomFactory_mod
+      use mapl_LocStreamGeomFactory_mod
       use mapl_XYGeomFactory_mod
       use mapl_EASEGeomFactory_mod
       class(GeomManager), intent(inout) :: this
 
       ! Load default factories
       type(LatLonGeomFactory) :: latlon_factory
       type(CubedSphereGeomFactory) :: cs_factory
+      type(LocStreamGeomFactory) :: locstream_factory
       type(XYGeomFactory) :: xy_factory
       type(EASEGeomFactory) :: ease_factory
 
       call this%add_factory(cs_factory)
+      call this%add_factory(locstream_factory)
       call this%add_factory(latlon_factory)
       call this%add_factory(xy_factory)
       call this%add_factory(ease_factory)