subroutine test_conserve_regrid_csgrid(isregular, rc)
#undef ESMF_METHOD
#define ESMF_METHOD "test_conserve_regrid_csgrid"
logical, intent(in) :: isregular
integer, intent(out) :: rc
logical :: correct
integer :: localrc
type(ESMF_Grid) :: srcGrid
type(ESMF_Grid) :: dstGrid
type(ESMF_Field) :: srcField
type(ESMF_Field) :: dstField
type(ESMF_Field) :: xdstField
type(ESMF_Field) :: errField
type(ESMF_Field) :: srcArea
type(ESMF_Field) :: dstArea
type(ESMF_Field) :: srcFracField
type(ESMF_Field) :: dstFracField
type(ESMF_Array) :: dstArray
type(ESMF_Array) :: errArray
type(ESMF_Array) :: srcArray
type(ESMF_RouteHandle) :: routeHandle
type(ESMF_VM) :: vm
real(ESMF_KIND_R8), pointer :: farrayPtrXC(:,:)
real(ESMF_KIND_R8), pointer :: farrayPtrYC(:,:)
real(ESMF_KIND_R8), pointer :: farrayPtr(:,:), farrayPtr2(:,:)
real(ESMF_KIND_R8), pointer :: xfarrayPtr(:,:)
real(ESMF_KIND_R8), pointer :: errfarrayPtr(:,:)
real(ESMF_KIND_R8), pointer :: srcAreaptr(:,:), dstAreaptr(:,:)
real(ESMF_KIND_R8), pointer :: srcFracptr(:,:), dstFracptr(:,:)
integer :: clbnd(2),cubnd(2)
integer :: fclbnd(2),fcubnd(2)
integer :: i1,i2,i3, index(2)
integer :: lDE, srclocalDECount, dstlocalDECount
real(ESMF_KIND_R8) :: coord(2),x,y,z
character(len=ESMF_MAXSTR) :: string
integer src_tile_size, dst_nx, dst_ny
real(ESMF_KIND_R8) :: lon, lat, theta, phi, relErr
real(ESMF_KIND_R8) :: coords(2), maxRelErr
integer :: localPet, petCount
real(ESMF_KIND_R8) :: dstIntegral, srcIntegral
real(ESMF_KIND_R8) :: localSrcIntegral(1), globalSrcIntegral(1)
real(ESMF_KIND_R8) :: localDstIntegral(1), globalDstIntegral(1)
real(ESMF_KIND_R8) :: localMaxRelErr(1), globalMaxRelErr(1)
real(ESMF_KIND_R8), parameter :: DEG2RAD = &
3.141592653589793_ESMF_KIND_R8/180.0_ESMF_KIND_R8
integer :: decomptile(2,6)
integer :: countsPerDEDim1(3,6), countsPerDEDim2(2,6)
integer :: i
! init success flag
correct=.true.
rc=ESMF_SUCCESS
! get pet info
call ESMF_VMGetGlobal(vm, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMGet(vm, petCount=petCount, localPet=localpet, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Establish the resolution of the grids
src_tile_size=20
dst_nx = 47
dst_ny = 47
if (isregular) then
! Set up decomposition for src Grid
decomptile(:,1)=(/2,2/)
decomptile(:,2)=(/2,2/)
decomptile(:,3)=(/2,2/)
decomptile(:,4)=(/2,2/)
decomptile(:,5)=(/2,2/)
decomptile(:,6)=(/2,2/)
! Create Src Grid
srcGrid=ESMF_GridCreateCubedSphere(tileSize=src_tile_size, &
regDecompPTile=decomptile, &
staggerLocList = (/ESMF_STAGGERLOC_CORNER, ESMF_STAGGERLOC_CENTER/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
else
do i=1,3
CountsPerDeDim1(:,i)=(/10,5,5/)
CountsPerDeDim1(:,i+3)=(/10,10,0/)
CountsPerDeDim2(:,i)=(/12,8/)
CountsPerDeDim2(:,i+3)=(/20,0/)
enddo
! Create Src Grid
srcGrid=ESMF_GridCreateCubedSphere(src_tile_size, &
CountsPerDeDim1, CountsPerDeDim2, &
staggerLocList = (/ESMF_STAGGERLOC_CORNER, ESMF_STAGGERLOC_CENTER/), &
indexflag = ESMF_INDEX_GLOBAL, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
endif
! create src fields
srcField = ESMF_FieldCreate(srcGrid, typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
srcFracField = ESMF_FieldCreate(srcGrid, typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
srcArea = ESMF_FieldCreate(srcGrid, typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="source", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! srcArray
call ESMF_FieldGet(srcField, array=srcArray, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get number of local DEs
call ESMF_GridGet(srcGrid, localDECount=srclocalDECount, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get memory and set fields for src
do lDE=0,srclocalDECount-1
! get src pointer
call ESMF_FieldGet(srcField, lDE, farrayPtr, &
computationalLBound=clbnd, computationalUBound=cubnd, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!! set coords, interpolated function
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
! Get coords
call ESMF_GridGetCoord(srcGrid, staggerloc=ESMF_STAGGERLOC_CENTER, &
localDE=lDE, index=(/i1,i2/), coord=coords, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! init exact answer
lon = coords(1)
lat = coords(2)
! Set the source to be a function of the x,y,z coordinate
theta = DEG2RAD*(lon)
phi = DEG2RAD*(90.-lat)
x = cos(theta)*sin(phi)
y = sin(theta)*sin(phi)
z = cos(phi)
! set src data
farrayPtr(i1,i2) = x+y+z+15.0
! This one seems to do a weird thing around the pole with a cubed sphere
! farrayPtr(i1,i2) = 2. + cos(theta)**2.*cos(2.*phi)
enddo
enddo
enddo ! lDE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Destination grid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Create Dst Grid
dstGrid=ESMF_GridCreate1PeriDimUfrm(maxIndex=(/dst_nx,dst_ny/), &
minCornerCoord=(/0.0_ESMF_KIND_R8,-90.0_ESMF_KIND_R8/), &
maxCornerCoord=(/360.0_ESMF_KIND_R8,90.0_ESMF_KIND_R8/), &
staggerLocList=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Create dst fields
dstField = ESMF_FieldCreate(dstGrid, typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="dest", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
xdstField = ESMF_FieldCreate(dstGrid, typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="xdest", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
errField = ESMF_FieldCreate(dstGrid, typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="xdest", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
dstFracField = ESMF_FieldCreate(dstGrid, typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="dest", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
dstArea = ESMF_FieldCreate(dstGrid, typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="dest", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get arrays
call ESMF_FieldGet(dstField, array=dstArray, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldGet(errField, array=errArray, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get number of local DEs
call ESMF_GridGet(dstGrid, localDECount=dstlocalDECount, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get memory and set fields for dst
do lDE=0,dstlocalDECount-1
! get dst pointer
call ESMF_FieldGet(dstField, lDE, farrayPtr, &
computationalLBound=clbnd, computationalUBound=cubnd, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get exact dst pointer
call ESMF_FieldGet(xdstField, lDE, xfarrayPtr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!! dst data
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
! Get coords
call ESMF_GridGetCoord(dstGrid, staggerloc=ESMF_STAGGERLOC_CENTER, &
localDE=lDE, index=(/i1,i2/), coord=coords, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! init exact answer
lon = coords(1)
lat = coords(2)
! Set the source to be a function of the x,y,z coordinate
theta = DEG2RAD*(lon)
phi = DEG2RAD*(90.-lat)
x = cos(theta)*sin(phi)
y = sin(theta)*sin(phi)
z = cos(phi)
! farrayPtr(i1,i2) = 2. + cos(theta)**2.*cos(2.*phi)
! set exact dst data
xfarrayPtr(i1,i2) = x+y+z+15.0
! This one seems to do a weird thing around the pole with a cubed sphere
!xfarrayPtr(i1,i2) = 2. + cos(theta)**2.*cos(2.*phi)
! initialize destination field
farrayPtr(i1,i2)=0.0
enddo
enddo
enddo ! lDE
#if 0
call ESMF_GridWriteVTK(srcGrid,staggerloc=ESMF_STAGGERLOC_CENTER, &
filename="srcGrid", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
#endif
!!! Regrid forward from the A grid to the B grid
! Regrid store
call ESMF_FieldRegridStore( &
srcField, &
dstField=dstField, &
routeHandle=routeHandle, &
unmappedAction=ESMF_UNMAPPEDACTION_ERROR, &
regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
dstFracField=dstFracField, &
srcFracField=srcFracField, &
normType=ESMF_NORMTYPE_FRACAREA, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Do regrid
call ESMF_FieldRegrid(srcField, dstField, routeHandle, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegridRelease(routeHandle, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get the src cell areas
call ESMF_FieldRegridGetArea(srcArea, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get the dst cell areas
call ESMF_FieldRegridGetArea(dstArea, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! XMRKX
! Check results
maxRelErr=0.0
dstIntegral=0.0
do lDE=0,dstlocalDECount-1
call ESMF_FieldGet(dstField, lDE, farrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldGet(xdstField, lDE, xfarrayPtr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get dst area Field
call ESMF_FieldGet(dstArea, lDE, dstAreaptr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get frac Field
call ESMF_FieldGet(dstFracField, lDE, dstFracptr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldGet(errField, lDE, errfarrayPtr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!! make sure we're not using any bad points
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
! Compute total integral
! NOTE: DO need to include dstFrac here, because the frac has been included in the weights
dstIntegral = dstIntegral + dstFracptr(i1,i2)*dstAreaptr(i1,i2)*farrayPtr(i1,i2)
! dstIntegral = dstIntegral + dstAreaptr(i1,i2)
! Compute relative error
if (xfarrayPtr(i1,i2) .ne. 0.0) then
relErr=abs((farrayPtr(i1,i2)-xfarrayPtr(i1,i2))/xfarrayPtr(i1,i2))
else
relErr=abs(farrayPtr(i1,i2)-xfarrayPtr(i1,i2))
endif
! if working everything should be close to exact answer
if (relErr .gt. 0.005) then
correct=.false.
!write(*,*) "relErr=",relErr,farrayPtr(i1,i2),xfarrayPtr(i1,i2)
endif
! Calc max
if (relErr > maxRelErr) then
maxRelErr=relErr
endif
! put in error field
errfarrayPtr(i1,i2)=relErr
enddo
enddo
enddo ! lDE
! Compute src integral
srcIntegral=0.0
do lDE=0,srcLocalDECount-1
! get src pointer
call ESMF_FieldGet(srcField, lDE, farrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get src Field
call ESMF_FieldGet(srcArea, lDE, srcAreaptr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get frac Field
call ESMF_FieldGet(srcFracField, lDE, srcFracptr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
srcIntegral = srcIntegral + srcFracptr(i1,i2)*srcAreaptr(i1,i2)*farrayPtr(i1,i2)
! srcIntegral = srcIntegral + srcAreaptr(i1,i2)
enddo
enddo
enddo ! lDE
#if 0
call ESMF_GridWriteVTK(srcGrid,staggerloc=ESMF_STAGGERLOC_CENTER, &
filename="srcGrid", array1=srcArray, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridWriteVTK(dstGrid,staggerloc=ESMF_STAGGERLOC_CENTER, &
filename="dstGrid", array1=dstArray, array2=errArray, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
#endif
! Get info across PETs
globalSrcIntegral(1) = 0.0
globalDstIntegral(1) = 0.0
localSrcIntegral(1)=srcIntegral
call ESMF_VMAllReduce(vm, localSrcIntegral, globalSrcIntegral, 1, &
ESMF_REDUCE_SUM, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
localDstIntegral(1)=dstIntegral
call ESMF_VMAllReduce(vm, localDstIntegral, globalDstIntegral, 1, &
ESMF_REDUCE_SUM, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
localMaxRelErr(1)=maxRelErr
call ESMF_VMAllReduce(vm, localMaxRelErr, globalMaxRelErr, 1, &
ESMF_REDUCE_MAX, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Uncomment these calls to see regrid results
#if 0
if (localPet == 0) then
write(*,*) "=== Cubed Sphere Grid ==="
write(*,*) "Conservation:"
write(*,*) "Rel Error = ", ABS(globalDstIntegral(1)-globalSrcIntegral(1))/globalSrcIntegral(1)
write(*,*) "SRC Integral = ", globalSrcIntegral(1)
write(*,*) "DST Integral = ", globalDstIntegral(1)
write(*,*) " "
write(*,*) "Interpolation:"
write(*,*) "Max Error = ", globalMaxRelErr(1)
write(*,*) " "
endif
#endif
! return answer based on values
if (ABS(globalDstIntegral(1)-globalSrcIntegral(1))/globalSrcIntegral(1) > 1.0E-14) correct=.false.
if (globalMaxRelErr(1) > 10E-2) correct=.false.
! Destroy the Fields
call ESMF_FieldDestroy(srcField, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldDestroy(dstField, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Free the grids
call ESMF_GridDestroy(srcGrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridDestroy(dstGrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! return answer based on correct flag
if (correct) then
rc=ESMF_SUCCESS
else
rc=ESMF_FAILURE
endif
end subroutine test_conserve_regrid_csgrid