Changes to return correct dimensions for the MPI grid from the OpenMP subgrids

CHANGELOG.md

@@ -10,10 +10,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 ### Changed
+
 - Updated ExtData so that if files are missing in a sequence the last value will be perisisted if one has not chosen `exact` option
 
 ### Fixed
 
+- Changes were made to add attributes to the subgrids (i.e. created by dividing the MPI subdomain into smaller subdomains equal to the number of OpenMP threads) such that the correct dimensions for the MPI subdomain could be retrieved from the subgrids where ever needed.
+
 ### Removed
 
 ### Deprecated

base/Base/Base_Base_implementation.F90

@@ -1557,11 +1557,24 @@ module subroutine MAPL_GRID_INTERIOR(GRID,I1,IN,J1,JN)
     integer                               :: gridRank
     integer,               allocatable    :: localDeToDeMap(:)
     integer :: rc
+    logical :: isPresent
+    integer :: global_grid_info(10)
 
     i1=-1
     j1=-1
     in=-1
     jn=-1
+
+    call ESMF_AttributeGet(grid, name="GLOBAL_GRID_INFO", isPresent=isPresent, _RC)
+    if (isPresent) then
+      call ESMF_AttributeGet(grid, name="GLOBAL_GRID_INFO", valueList=global_grid_info, _RC)
+      I1 = global_grid_info(7)
+      IN = global_grid_info(8)
+      j1 = global_grid_info(9)
+      JN = global_grid_info(10)
+      _RETURN(_SUCCESS)
+    end if
+
     call ESMF_GridGet    (GRID, dimCount=gridRank, distGrid=distGrid, _RC)
     call ESMF_DistGridGet(distGRID, delayout=layout, _RC)
     call ESMF_DELayoutGet(layout, deCount = nDEs, localDeCount=localDeCount,_RC)
@@ -2138,6 +2151,24 @@ module subroutine MAPL_GridGetInterior(GRID,I1,IN,J1,JN)
     integer                               :: deId
     integer                               :: gridRank
     integer                               :: rc
+    logical                               :: isPresent
+    integer                               :: global_grid_info(10)
+
+    i1=-1
+    j1=-1
+    in=-1
+    jn=-1
+
+    call ESMF_AttributeGet(grid, name="GLOBAL_GRID_INFO", isPresent=isPresent, _RC)
+    if (isPresent) then
+      call ESMF_AttributeGet(grid, name="GLOBAL_GRID_INFO", valueList=global_grid_info, _RC)
+      I1 = global_grid_info(7)
+      IN = global_grid_info(8)
+      j1 = global_grid_info(9)
+      JN = global_grid_info(10)
+      _RETURN(_SUCCESS)
+    end if
+
 
     call ESMF_GridGet    (GRID, dimCount=gridRank, distGrid=distGrid, _RC)
     call ESMF_DistGridGet(distGRID, delayout=layout, _RC)
@@ -2627,9 +2658,7 @@ module subroutine MAPL_GetHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, rc)
        tmp_lats = latR8
     end if
 
-!AOO change tusing GridType atribute    if (im_world*6==jm_world) then
-    call ESMF_AttributeGet(grid, name='GridType', value=grid_type, _RC)
-    if(trim(grid_type) == "Cubed-Sphere") then
+    if (im_world*6==jm_world) then
 
       call MAPL_GetGlobalHorzIJIndex(npts, II, JJ, lon=lon, lat=lat, lonR8=lonR8, latR8=latR8, Grid=Grid, _RC)
 
@@ -2868,28 +2897,33 @@ function grid_is_ok(grid) result(OK)
        type(ESMF_Grid), intent(inout) :: grid
        logical :: OK
        integer :: I1, I2, J1, J2, j
-       real(ESMF_KIND_R8), pointer :: corner_lons(:,:), corner_lats(:,:)
+       real(ESMF_KIND_R8), allocatable :: corner_lons(:,:), corner_lats(:,:)
        real(ESMF_KIND_R8), allocatable :: lonRe(:), latRe(:)
        real(ESMF_KIND_R8), allocatable :: accurate_lat(:), accurate_lon(:)
        real(ESMF_KIND_R8) :: stretch_factor, target_lon, target_lat, shift0
        real :: tolerance
+       integer :: local_dims(3)
 
        tolerance = epsilon(1.0)
        call MAPL_GridGetInterior(grid,I1,I2,J1,J2)
+       call MAPL_GridGet(grid, localCellCountPerDim=local_dims, _RC)
        OK = .true.
        ! check the edge of face 1 along longitude
-       call ESMF_GridGetCoord(grid,localDE=0,coordDim=1,staggerloc=ESMF_STAGGERLOC_CORNER, &
-            farrayPtr=corner_lons, rc=status)
-       call ESMF_GridGetCoord(grid,localDE=0,coordDim=2,staggerloc=ESMF_STAGGERLOC_CORNER, &
-            farrayPtr=corner_lats, rc=status)
+       !call ESMF_GridGetCoord(grid,localDE=0,coordDim=1,staggerloc=ESMF_STAGGERLOC_CORNER, &
+       !     farrayPtr=corner_lons, _RC)
+       !call ESMF_GridGetCoord(grid,localDE=0,coordDim=2,staggerloc=ESMF_STAGGERLOC_CORNER, &
+       !     farrayPtr=corner_lats, _RC)
+       allocate(corner_lons(local_dims(1)+1, local_dims(2)+1))
+       allocate(corner_lats(local_dims(1)+1, local_dims(2)+1))
+       call MAPL_GridGetCorners(grid, corner_lons, corner_lats, _RC)
 
 
        if ( I1 == 1 .and. J1 == 1 ) then
-          allocate(lonRe(j2-j1+1),      latRe(j2-j1+1))
-          call MAPL_Reverse_Schmidt(grid, stretched, J2-J1+1, lonR8=corner_lons(1,:), &
-                                   latR8=corner_lats(1,:), lonRe=lonRe, latRe=latRe, _RC)  
+          allocate(lonRe(local_dims(2)),      latRe(local_dims(2)))
+          call MAPL_Reverse_Schmidt(grid, stretched, local_dims(2), lonR8=corner_lons(1,1:local_dims(2)), &
+                                   latR8=corner_lats(1,1:local_dims(2)), lonRe=lonRe, latRe=latRe, _RC)  
 
-          allocate(accurate_lon(j2-j1+1), accurate_lat(j2-j1+1))
+          allocate(accurate_lon(local_dims(2)), accurate_lat(local_dims(2)))
 
           shift0 = shift
           if (stretched) shift0 = 0

base/MaplGrid.F90

@@ -265,10 +265,19 @@ subroutine MAPL_GridGet(GRID, globalCellCountPerDim, localCellCountPerDim, layou
       type(ESMF_DistGrid) :: distGrid
       integer, allocatable :: maxindex(:,:),minindex(:,:)
       integer, pointer :: ims(:),jms(:)
+      integer :: global_grid_info(10)
 
       pglobal = present(globalCellCountPerDim)
       plocal  = present(localCellCountPerDim)
 
+      call ESMF_AttributeGet(grid, name="GLOBAL_GRID_INFO", isPresent=isPresent, _RC)
+      if (isPresent) then
+        call ESMF_AttributeGet(grid, name="GLOBAL_GRID_INFO", valueList=global_grid_info, _RC)
+        if (pglobal) globalCellCountPerDim = global_grid_info(1:3)
+        if (plocal)  localCellCountPerDim = global_grid_info(4:6)
+        _RETURN(_SUCCESS)
+      end if
+
       if (pglobal .or. plocal) then
          call ESMF_GridGet(grid, dimCount=gridRank, _RC)
 

generic/OpenMP_Support.F90

@@ -6,6 +6,7 @@ module MAPL_OpenMP_Support
     use MAPL_maplgrid
     use MAPL_ExceptionHandling
     use mapl_KeywordEnforcerMod
+    use MAPL_BaseMod, only : MAPL_Grid_Interior
     !$ use omp_lib
 
     implicit none
@@ -83,7 +84,7 @@ function make_subgrids_from_bounds(primary_grid, bounds, unusable, rc) result(su
         type(Interval), intent(in) :: bounds(:)
         class(KeywordEnforcer), optional, intent(in) :: unusable
         integer, optional, intent(out) :: rc
-        integer :: local_count(3)
+        integer :: local_count(3), global_count(3)
         integer :: status
         integer :: petMap(1,1,1)
         integer :: myPet, section, i, j, k, count, size_
@@ -102,7 +103,7 @@ function make_subgrids_from_bounds(primary_grid, bounds, unusable, rc) result(su
         !end do
 
         allocate(subgrids(size(bounds)))
-        call MAPL_GridGet(primary_grid,localcellcountPerDim=local_count, _RC)
+        call MAPL_GridGet(primary_grid,localcellcountPerDim=local_count, globalCellCountPerDim=global_count, _RC)
         call ESMF_VMGetCurrent(vm, _RC)
         call ESMF_VMGet(vm, localPET=myPET, _RC)
 
@@ -175,6 +176,22 @@ function make_subgrids_from_bounds(primary_grid, bounds, unusable, rc) result(su
                 itemCount = count, valueList=lons1d, _RC)
            call ESMF_AttributeSet(subgrids(i), name='GridCornerLats:', &
                 itemCount = count, valueList=lats1d, _RC)
+           block
+              integer :: global_grid_info(10)
+              integer :: i1,i2,j1,j2
+              call MAPL_Grid_Interior(primary_grid,i1,i2,j1,j2)
+              global_grid_info(1:3) = global_count
+              !global_grid_info(4:6) = local_count
+              global_grid_info(4) = size(new_lons,1)
+              global_grid_info(5) = size(new_lons,2)
+              global_grid_info(6) = local_count(3)
+              global_grid_info(7) = i1
+              global_grid_info(8) = i2
+              global_grid_info(9) = j1 + bounds(i)%min - 1
+              global_grid_info(10) = j1 + bounds(i)%max - 1
+              call ESMF_AttributeSet(subgrids(i), name="GLOBAL_GRID_INFO",  &
+                   itemCount=10, valueList=global_grid_info, _RC)
+           end block
 
             deallocate(lons1d, lats1d)
             deallocate(new_corner_lons, new_corner_lats)