subroutine MAPL_GridCoordAdjust(GRID, LOCSTREAM, RC)
type(ESMF_Grid), intent(INout ) :: Grid
type(MAPL_LocStream), intent(IN ) :: Locstream
integer, optional, intent(OUT) :: RC
! local vars
!------------
integer :: STATUS
integer :: NGRIDS
integer :: I, J, N
integer :: IM, JM
logical :: found
integer :: COUNTS(3)
integer :: NT, IG
character(len=MAPL_TileNameLength) :: GRIDNAME
character(len=MAPL_TileNameLength), pointer :: GNAMES(:)
real(ESMF_KIND_R8) :: X, Y, W
real(ESMF_KIND_R8), allocatable :: sumw(:,:), sumxw(:,:), sumyw(:,:)
real(ESMF_KIND_R8), pointer :: gridx(:,:), gridy(:,:)
! get grid name
call ESMF_GridGet(grid, name=gridname, rc=status)
_VERIFY(STATUS)
call MAPL_LocstreamGet(LOCSTREAM, GRIDNAMES = GNAMES, RC=STATUS)
_VERIFY(STATUS)
! query loc_in for ngrids
ngrids = size(gnames)
_ASSERT(ngrids==2,'needs informative message')
! validate that gridname_in is there
found = .false.
DO I = 1, NGRIDS
IF (GNAMES(I) == GRIDNAME) THEN
FOUND = .TRUE.
exit
ENDIF
ENDDO
_ASSERT(FOUND,'needs informative message')
! get id of the grid we just found
IG = I
_ASSERT(IG == LocStream%Ptr%Current_Tiling,'needs informative message')
! get IM, JM and IM_WORLD, JM_WORLD
call MAPL_GridGet(GRID, localCellCountPerDim=COUNTS, RC=STATUS)
_VERIFY(STATUS)
IM = COUNTS(1)
JM = COUNTS(2)
! Retrieve the coordinates so we can set them
call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, &
staggerloc=ESMF_STAGGERLOC_CENTER, &
farrayPtr=gridX, rc=status)
_VERIFY(STATUS)
call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, &
staggerloc=ESMF_STAGGERLOC_CENTER, &
farrayPtr=gridY, rc=status)
_VERIFY(STATUS)
allocate(sumxw(IM, JM), sumyw(IM, JM), sumw (IM, JM), stat=status)
_VERIFY(STATUS)
SUMW = 0.0
SUMXW = 0.0
SUMYW = 0.0
NT = LOCSTREAM%Ptr%NT_Local
! loop over tiles
DO N = 1, NT
I = LOCSTREAM%Ptr%Local_IndexLocation(N)%I
J = LOCSTREAM%Ptr%Local_IndexLocation(N)%J
W = LOCSTREAM%Ptr%Local_IndexLocation(N)%W
X = locstream%Ptr%Local_GeoLocation(N)%X
Y = locstream%Ptr%Local_GeoLocation(N)%Y
SUMW(I,J) = SUMW(I,J) + W
SUMXW(I,J) = SUMXW(I,J) + X * W
SUMYW(I,J) = SUMYW(I,J) + Y * W
END DO
WHERE (SUMW == 0.0)
SUMXW = MAPL_UNDEF
SUMYW = MAPL_UNDEF
ELSEWHERE
SUMXW = SUMXW / SUMW
SUMYW = SUMYW / SUMW
! Make sure the longitudes are between -180 and 180 degrees
SUMXW = mod(SUMXW + 72180._REAL64,360._REAL64) - 180._REAL64 ! -180<= lon0 <180
! Convert to radians
SUMXW = SUMXW * (MAPL_PI_R8)/180._REAL64
SUMYW = SUMYW * (MAPL_PI_R8)/180._REAL64
END WHERE
! Modify grid coordinates
!------------------------
GRIDX = SUMXW
GRIDY = SUMYW
! Clean-up
!---------
deallocate(sumw)
deallocate(sumyw)
deallocate(sumxw)
! All done
!---------
_RETURN(ESMF_SUCCESS)
end subroutine MAPL_GridCoordAdjust