subroutine test_MeshToMesh_2nd(rc)
#undef ESMF_METHOD
#define ESMF_METHOD "test_MeshToMesh_2nd"
integer, intent(out) :: rc
logical :: itrp
logical :: csrv
integer :: localrc
type(ESMF_Mesh) :: srcMesh
type(ESMF_Mesh) :: dstMesh
type(ESMF_XGrid) :: xgrid
type(ESMF_Field) :: srcField
type(ESMF_Field) :: dstField
type(ESMF_Field) :: xdstField
type(ESMF_Field) :: xField
type(ESMF_Field) :: srcAreaField, dstAreaField
type(ESMF_Field) :: srcFracField, dstFracField
type(ESMF_RouteHandle) :: StoXrouteHandle
type(ESMF_RouteHandle) :: XtoDrouteHandle
type(ESMF_RouteHandle) :: routeHandle
type(ESMF_ArraySpec) :: arrayspec
type(ESMF_VM) :: vm
real(ESMF_KIND_R8), pointer :: srcFarrayPtr(:), dstFarrayPtr(:), xdstFarrayPtr(:)
real(ESMF_KIND_R8), pointer :: srcAreaPtr(:), dstAreaPtr(:)
real(ESMF_KIND_R8), pointer :: srcFracPtr(:), dstFracPtr(:)
integer :: clbnd(1),cubnd(1)
integer :: i1,i2,i3
real(ESMF_KIND_R8) :: x,y,z
real(ESMF_KIND_R8) :: lat, lon, phi, theta
real(ESMF_KIND_R8),parameter :: &
DEG2RAD = 3.141592653589793_ESMF_KIND_R8/180.0_ESMF_KIND_R8
integer :: localPet, petCount
real(ESMF_KIND_R8) :: srcmass(1), dstmass(1), srcmassg(1), dstmassg(1)
real(ESMF_KIND_R8) :: maxerror(1), minerror(1), error
real(ESMF_KIND_R8) :: maxerrorg(1), minerrorg(1), errorg
real(ESMF_KIND_R8) :: errorTot, errorTotG, dstVal
integer :: num_errorTot
real(ESMF_KIND_R8) :: l_errorTot(1),g_errorTot(1)
integer :: l_num_errorTot(1), g_num_errorTot(1)
integer :: numOwnedElems
real(ESMF_KIND_R8), pointer :: ownedElemCoords(:)
! result code
integer :: finalrc
! Init to success
rc=ESMF_SUCCESS
itrp=.true.
csrv=.true.
! 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
! If we don't have 1 or 2 PETS then exit unsuccessfully
if ((petCount .ne. 1) .and. (petCount .ne. 2)) then
rc=ESMF_FAILURE
return
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!! Setup Source !!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Create Source Mesh
call CreateTestMesh2x2_1(srcMesh, rc=localrc) ! Non-easy element create
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Array spec for fields
call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Create source field
srcField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Create source area field
srcAreaField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source_area", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Create source frac field
srcFracField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source_frac", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Load test data into the source Field
! Should only be 1 localDE
call ESMF_FieldGet(srcField, 0, srcFarrayPtr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Set interpolated function
call ESMF_MeshGet(srcMesh, numOwnedElements=numOwnedElems, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Allocate space for coordinates
allocate(ownedElemCoords(2*numOwnedElems))
! Set interpolated function
call ESMF_MeshGet(srcMesh, ownedElemCoords=ownedElemCoords, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! loop through and set field
do i1=1,numOwnedElems
! Get coords
lon=ownedElemCoords(2*i1-1)
lat=ownedElemCoords(2*i1)
! Set source function
theta = DEG2RAD*(lon)
phi = DEG2RAD*(90.-lat)
x = cos(theta)*sin(phi)
y = sin(theta)*sin(phi)
z = cos(phi)
srcFarrayPtr(i1) = x+y+z
!srcFarrayPtr(i1) = 1.0
enddo
! Deallocate space for coordinates
deallocate(ownedElemCoords)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!! Setup Destination !!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Create Destination Mesh
call CreateTestMesh2x2_2(dstMesh, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Array spec
call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Create dest. field
dstField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Create dest. area field
dstAreaField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest_area", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Create dest. frac field
dstFracField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest_frac", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Create exact dest. field
xdstField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="xdest", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Init destination field to 0.0
! Should only be 1 localDE
call ESMF_FieldGet(dstField, 0, dstFarrayPtr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Init exact destination field
! Should only be 1 localDE
call ESMF_FieldGet(xdstField, 0, xdstFarrayPtr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Set number of points in destination mesh
call ESMF_MeshGet(dstMesh, numOwnedElements=numOwnedElems, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Allocate space for coordinates
allocate(ownedElemCoords(2*numOwnedElems))
! Set exact destination field
call ESMF_MeshGet(dstMesh, ownedElemCoords=ownedElemCoords, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! loop through and set xfield
do i1=1,numOwnedElems
! Get coords
lon=ownedElemCoords(2*i1-1)
lat=ownedElemCoords(2*i1)
! Set exact dest function
theta = DEG2RAD*(lon)
phi = DEG2RAD*(90.-lat)
x = cos(theta)*sin(phi)
y = sin(theta)*sin(phi)
z = cos(phi)
xdstFarrayPtr(i1) = x+y+z
! xdstFarrayPtr(i1) = 1.0
! Init destination field to 0.0
dstFarrayPtr(i1)=0.0
enddo
! Deallocate space for coordinates
deallocate(ownedElemCoords)
#if 0
call ESMF_MeshWrite(srcMesh,"srcMesh")
call ESMF_MeshWrite(dstMesh,"dstMesh")
#endif
#define USE_XGRID
#ifdef USE_XGRID
!write(*,*) "Using XGrid"
! Create XGrid
xgrid = ESMF_XGridCreate(sideAMesh=(/srcMesh/), &
sideBMesh=(/dstMesh/), &
storeOverlay = .true., &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Field on XGrid
xField = ESMF_FieldCreate(xgrid, arrayspec, &
name="xfield", rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Regrid store
call ESMF_FieldRegridStore( &
xgrid, &
srcField, &
dstField=xField, &
routeHandle=StoXrouteHandle, &
regridmethod=ESMF_REGRIDMETHOD_CONSERVE_2ND, &
srcFracField=srcFracField, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegridStore( &
xgrid, &
xField, &
dstField=dstField, &
routeHandle=XToDrouteHandle, &
regridmethod=ESMF_REGRIDMETHOD_CONSERVE_2ND, &
dstFracField=dstFracField, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Do regrid
call ESMF_FieldRegrid(srcField, xField, StoXrouteHandle, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegrid(xField, dstField, XToDrouteHandle, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Release routehandles
call ESMF_FieldRegridrelease(StoXrouteHandle, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegridRelease(XtoDrouteHandle, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
#else
!write(*,*) "NOT Using XGrid"
!!! Regrid forward from the A grid to the B grid
! Regrid store
call ESMF_FieldRegridStore( &
srcField, &
dstField=dstField, &
routeHandle=routeHandle, &
regridmethod=ESMF_REGRIDMETHOD_CONSERVE_2ND, &
dstFracField=dstFracField, &
srcFracField=srcFracField, &
unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, &
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
#endif
! Get the integration weights
call ESMF_FieldRegridGetArea(srcAreaField, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get the integration weights
call ESMF_FieldRegridGetArea(dstAreaField, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Check if the values are close
minerror(1) = 100000.
maxerror(1) = 0.
error = 0.
errorTot=0.0
num_errorTot=0
dstmass = 0.
! get dst Field
call ESMF_FieldGet(dstField, 0, dstFarrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get exact destination Field
call ESMF_FieldGet(xdstField, 0, xdstFarrayPtr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get dst area Field
call ESMF_FieldGet(dstAreaField, 0, dstAreaPtr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get frac Field
call ESMF_FieldGet(dstFracField, 0, dstFracptr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! destination grid
!! check relative error
do i1=clbnd(1),cubnd(1)
! This is WRONG, shouldn't include Frac
! dstmass = dstmass + dstFracptr(i1,i2)*dstAreaptr(i1)*fptr(i1)
! Instead do this
dstmass = dstmass + dstAreaptr(i1)*dstFarrayPtr(i1)
! If this destination cell isn't covered by a sig. amount of source, then don't compute error on it.
if (dstFracPtr(i1) .lt. 0.1) cycle
! write(*,*) i1,"::",dstFarrayPtr(i1),xdstFarrayPtr(i1)
! Since fraction isn't included in weights in this case, use it to modify dstField value, so
! that it's correct for partial cells
if (dstFracPtr(i1) .ne. 0.0) then
dstVal=dstFarrayPtr(i1)/dstFracptr(i1)
else
dstVal=dstFarrayPtr(i1)
endif
if (xdstFarrayPtr(i1) .ne. 0.0) then
error=ABS(dstVal - xdstFarrayPtr(i1))/ABS(xdstFarrayPtr(i1))
else
error=ABS(dstVal - xdstFarrayPtr(i1))
endif
! total error
errorTot=errorTot+error
num_errorTot=num_errorTot+1
! min max error
if (error > maxerror(1)) then
maxerror(1) = error
endif
if (error < minerror(1)) then
minerror(1) = error
endif
enddo
srcmass(1) = 0.
! get src pointer
call ESMF_FieldGet(srcField, 0, srcFarrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get src Field
call ESMF_FieldGet(srcAreaField, 0, srcAreaptr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! get frac Field
call ESMF_FieldGet(srcFracField, 0, srcFracptr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i1=clbnd(1),cubnd(1)
srcmass(1) = srcmass(1) + srcFracptr(i1)*srcAreaptr(i1)*srcFarrayPtr(i1)
enddo
! Init integrals
srcmassg(1) = 0.
dstmassg(1) = 0.
call ESMF_VMAllReduce(vm, srcmass, srcmassg, 1, ESMF_REDUCE_SUM, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMAllReduce(vm, dstmass, dstmassg, 1, ESMF_REDUCE_SUM, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMAllReduce(vm, maxerror, maxerrorg, 1, ESMF_REDUCE_MAX, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMAllReduce(vm, minerror, minerrorg, 1, ESMF_REDUCE_MIN, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
l_errorTot(1)=errorTot
call ESMF_VMAllReduce(vm, l_errorTot, g_errorTot, 1, ESMF_REDUCE_SUM, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
l_num_errorTot(1)=num_errorTot
call ESMF_VMAllReduce(vm, l_num_errorTot, g_num_errorTot, 1, ESMF_REDUCE_SUM, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! return answer based on correct flags
csrv = .false.
if (ABS(dstmassg(1)-srcmassg(1))/srcmassg(1) < 1.0E-14) csrv = .true.
itrp = .false.
if (maxerrorg(1) < 1.5E-2) itrp = .true.
! Uncomment these calls to see some actual regrid results
if (localPet == 0) then
write(*,*)
write(*,*) "=== Second Order Conservative Mesh to Mesh via XGrid ==="
write(*,*) "Conservation:"
write(*,*) "Rel Error = ", ABS(dstmassg(1)-srcmassg(1))/srcmassg(1)
write(*,*) "SRC mass = ", srcmassg(1)
write(*,*) "DST mass = ", dstmassg(1)
write(*,*) " "
write(*,*) "Interpolation:"
write(*,*) "Max Error = ", maxerrorg(1)
write(*,*) "Min Error = ", minerrorg(1)
write(*,*) "Avg Error = ", g_errorTot(1)/g_num_errorTot(1)
write(*,*) " "
endif
! 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
call ESMF_FieldDestroy(srcAreaField, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldDestroy(dstAreaField, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldDestroy(srcFracField, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldDestroy(dstFracField, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldDestroy(xdstField, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Free the meshes
call ESMF_MeshDestroy(srcMesh, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_MeshDestroy(dstMesh, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! If either interpolation or conservation errors
! are too large, then return failure
if (.not. itrp .or. .not. csrv) then
rc=ESMF_FAILURE
endif
end subroutine test_MeshToMesh_2nd