subroutine test_RegridCsrv3DCartMesh(itrp, csrv, rc)
logical, intent(out) :: itrp
logical, intent(out) :: csrv
integer, intent(out) :: rc
integer :: localrc
type(ESMF_Mesh) :: srcMesh
type(ESMF_Mesh) :: dstMesh
type(ESMF_Field) :: srcField
type(ESMF_Field) :: dstField
type(ESMF_Field) :: xdstField
type(ESMF_Field) :: srcAreaField, dstAreaField
type(ESMF_Field) :: srcFracField, dstFracField
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
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
integer, pointer :: nodeIds(:),nodeOwners(:)
real(ESMF_KIND_R8), pointer :: nodeCoords(:)
integer, pointer :: elemIds(:),elemTypes(:),elemConn(:),elemMask(:)
integer :: numNodes, numElems
integer :: iconn,inode
! result code
integer :: finalrc
! Init to success
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
! If we don't have 1 or 4 PETS then exit successfully
if ((petCount .ne. 1) .and. (petCount .ne. 4)) then
print*,'ERROR: test must be run using exactly 1 or 4 PETS - detected ',petCount
rc=ESMF_FAILURE
return
endif
! Setup source mesh
if (petCount .eq. 1) then
! Set number of nodes
numNodes=18
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/1,2,3,4,5,6,7,8,9, &
10,11,12,13,14,15,16,17,18/)
! Allocate and fill node coordinate array.
! Since this is a 3D Mesh the size is 3x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/0.0,0.0,0.0, & ! node id 1
1.0,0.0,0.0, & ! node id 2
2.0,0.0,0.0, & ! node id 3
0.0,1.0,0.0, & ! node id 4
1.0,1.0,0.0, & ! node id 5
2.0,1.0,0.0, & ! node id 6
0.0,2.0,0.0, & ! node id 7
1.0,2.0,0.0, & ! node id 8
2.0,2.0,0.0, & ! node id 9
0.0,0.0,1.0, & ! node id 10
1.0,0.0,1.0, & ! node id 11
2.0,0.0,1.0, & ! node id 12
0.0,1.0,1.0, & ! node id 13
1.0,1.0,1.0, & ! node id 14
2.0,1.0,1.0, & ! node id 15
0.0,2.0,1.0, & ! node id 16
1.0,2.0,1.0, & ! node id 17
2.0,2.0,1.0 /) ! node id 18
! Allocate and fill the node owner array.
! Since this Mesh is all on PET 0, it's just set to all 0.
allocate(nodeOwners(numNodes))
nodeOwners=0 ! everything on PET 0
! Set the number of each type of element, plus the total number.
numElems=4
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/1,2,3,4/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=ESMF_MESHELEMTYPE_HEX
! Allocate and fill the element mask array
allocate(elemMask(numElems))
elemMask=(/0,0,0,1/)
! Allocate and fill the element connection type array.
! Note that entries in this array refer to the
! positions in the nodeIds, etc. arrays and that
! the order and number of entries for each element
! reflects that given in the Mesh options
! section for the corresponding entry
! in the elemTypes array.
allocate(elemConn(8*numElems))
elemConn=(/1,2,5,4,10,11,14,13, & ! elem id 1
2,3,6,5,11,12,15,14, & ! elem id 2
4,5,8,7,13,14,17,16, & ! elem id 3
5,6,9,8,14,15,18,17/) ! elem id 4
else if (petCount .eq. 4) then
! Setup mesh data depending on PET
if (localPET .eq. 0) then !!! This part only for PET 0
! Set number of nodes
numNodes=8
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/1,2,4,5,10,11,13,14/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/0.0,0.0,0.0, & ! node id 1
1.0,0.0,0.0, & ! node id 2
0.0,1.0,0.0, & ! node id 4
1.0,1.0,0.0, & ! node id 5
0.0,0.0,1.0, & ! node id 10
1.0,0.0,1.0, & ! node id 11
0.0,1.0,1.0, & ! node id 13
1.0,1.0,1.0 /) ! node id 14
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 1
0, & ! node id 2
0, & ! node id 4
0, & ! node id 5
0, & ! node id 10
0, & ! node id 11
0, & ! node id 13
0/) ! node id 14
! Set the number of each type of element, plus the total number.
numElems=1
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/1/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_HEX/) ! elem id 1
! Allocate and fill the element mask array
allocate(elemMask(numElems))
elemMask=(/0/)
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(8*numElems))
elemConn=(/1,2,4,3,5,6,8,7/) ! elem id 1
else if (localPET .eq. 1) then !!! This part only for PET 1
! Set number of nodes
numNodes=8
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/2,3,5,6,11,12,14,15/)
! Allocate and fill node coordinate array.
! Since this is a 3D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/1.0,0.0,0.0, & ! node id 2
2.0,0.0,0.0, & ! node id 3
1.0,1.0,0.0, & ! node id 5
2.0,1.0,0.0, & ! node id 6
1.0,0.0,1.0, & ! node id 11
2.0,0.0,1.0, & ! node id 12
1.0,1.0,1.0, & ! node id 14
2.0,1.0,1.0/) ! node id 15
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 2
1, & ! node id 3
0, & ! node id 5
1, & ! node id 6
0, & ! node id 11
1, & ! node id 12
0, & ! node id 14
1/) ! node id 15
! Set the number of each type of element, plus the total number.
numElems=1
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/2/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_HEX/) ! elem id 1
! Allocate and fill the element mask array
allocate(elemMask(numElems))
elemMask=(/0/)
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(8*numElems))
elemConn=(/1,2,4,3,5,6,8,7/) ! elem id 1
else if (localPET .eq. 2) then !!! This part only for PET 2
! Set number of nodes
numNodes=8
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/4,5,7,8,13,14,16,17/)
! Allocate and fill node coordinate array.
! Since this is a 3D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/0.0,1.0,0.0, & ! node id 4
1.0,1.0,0.0, & ! node id 5
0.0,2.0,0.0, & ! node id 7
1.0,2.0,0.0, & ! node id 8
0.0,1.0,1.0, & ! node id 13
1.0,1.0,1.0, & ! node id 14
0.0,2.0,1.0, & ! node id 16
1.0,2.0,1.0/) ! node id 17
! Allocate and fill the node owner array.
! Since this Mesh is all on PET 0, it's just set to all 0.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 4
0, & ! node id 5
2, & ! node id 7
2, & ! node id 8
0, & ! node id 13
0, & ! node id 14
2, & ! node id 16
2/) ! node id 17
! Set the number of each type of element, plus the total number.
numElems=1
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/3/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_HEX/) ! elem id 1
! Allocate and fill the element mask array
allocate(elemMask(numElems))
elemMask=(/0/)
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(8*numElems))
elemConn=(/1,2,4,3,5,6,8,7/) ! elem id 1
else if (localPET .eq. 3) then !!! This part only for PET 3
! Set number of nodes
numNodes=8
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/5,6,8,9,14,15,17,18/)
! Allocate and fill node coordinate array.
! Since this is a 3D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/1.0,1.0,0.0, & ! node id 5
2.0,1.0,0.0, & ! node id 6
1.0,2.0,0.0, & ! node id 8
2.0,2.0,0.0, & ! node id 9
1.0,1.0,1.0, & ! node id 14
2.0,1.0,1.0, & ! node id 15
1.0,2.0,1.0, & ! node id 17
2.0,2.0,1.0/) ! node id 18
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 5
1, & ! node id 6
2, & ! node id 8
3, & ! node id 9
0, & ! node id 14
1, & ! node id 15
2, & ! node id 17
3/) ! node id 18
! Set the number of each type of element, plus the total number.
numElems=1
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/4/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_HEX/) ! elem id 1
! Allocate and fill the element mask array
allocate(elemMask(numElems))
elemMask=(/1/)
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(8*numElems))
elemConn=(/1,2,4,3,5,6,8,7/) ! elem id 1
endif
endif
! Create Mesh structure in 1 step
srcMesh=ESMF_MeshCreate(parametricDim=3,spatialDim=3, &
coordSys=ESMF_COORDSYS_CART, &
nodeIds=nodeIds, nodeCoords=nodeCoords, &
nodeOwners=nodeOwners, elementIds=elemIds,&
elementTypes=elemTypes, elementConn=elemConn, &
elementMask=elemMask, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Array spec for fields
call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=rc)
! Create source field
srcField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source area field
srcAreaField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source_area", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source frac field
srcFracField = ESMF_FieldCreate(srcMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="source_frac", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Load test data into the source Field
! Should only be 1 localDE
call ESMF_FieldGet(srcField, 0, srcFarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! set interpolated function
iconn=1
do i1=1,numElems
! Loop through nodes in elem
! to compute point in center
x=0.0
y=0.0
z=0.0
do i3=1,8 ! 8 is the number of nodes in each element
inode=3*(elemConn(iconn)-1)
x=x+nodeCoords(inode+1)
y=y+nodeCoords(inode+2)
z=z+nodeCoords(inode+3)
iconn=iconn+1
enddo
x=0.125*x
y=0.125*y
z=0.125*z
! Set source function
! (Set huge value for masked values, so it can be detected)
if (elemMask(i1) .eq. 0) then
srcFarrayPtr(i1) = 20.0+x+y+z
else
srcFarrayPtr(i1) = 100000000.0
endif
enddo
! For now, Easy set interpolated function
!srcFarrayPtr=1.0
! deallocate node data
deallocate(nodeIds)
deallocate(nodeCoords)
deallocate(nodeOwners)
! deallocate elem data
deallocate(elemIds)
deallocate(elemTypes)
deallocate(elemConn)
deallocate(elemMask)
!!!!!!!!!!!!!!! setup destination Mesh !!!!!!!!!!!!!!!!!
! setup source Mesh
if (petCount .eq. 1) then
! Set number of nodes
numNodes=10
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/1,2,3,4,5,6,7,8,9, &
10/)
! Allocate and fill node coordinate array.
! Since this is a 3D Mesh the size is 3x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/0.0,0.0,0.0, & ! node id 1
1.0,0.0,0.0, & ! node id 2
2.0,0.0,0.0, & ! node id 3
0.5,1.0,0.0, & ! node id 4
1.5,1.0,0.0, & ! node id 5
1.0,2.0,0.0, & ! node id 6
0.5,0.5,1.0, & ! node id 7
1.0,0.5,1.0, & ! node id 8
1.5,0.5,1.0, & ! node id 9
1.0,1.5,1.0/) ! node id 10
! Allocate and fill the node owner array.
! Since this Mesh is all on PET 0, it's just set to all 0.
allocate(nodeOwners(numNodes))
nodeOwners=0 ! everything on PET 0
! Set the number of each type of element, plus the total number.
numElems=4
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/1,2,3,4/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=ESMF_MESHELEMTYPE_TETRA
! Allocate and fill the element connection type array.
! Note that entries in this array refer to the
! positions in the nodeIds, etc. arrays and that
! the order and number of entries for each element
! reflects that given in the Mesh options
! section for the corresponding entry
! in the elemTypes array.
allocate(elemConn(4*numElems))
elemConn=(/1,2,7,4, & ! elem id 1
2,3,9,5, & ! elem id 2
2,5,8,4, & ! elem id 3
4,5,10,6/) ! elem id 4
else if (petCount .eq. 4) then
! Setup mesh data depending on PET
if (localPET .eq. 0) then !!! This part only for PET 0
! Set number of nodes
numNodes=4
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/1,2,4,7/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/0.0,0.0,0.0, & ! node id 1
1.0,0.0,0.0, & ! node id 2
0.5,1.0,0.0, & ! node id 4
0.5,0.5,1.0/) ! node id 7
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 1
0, & ! node id 2
0, & ! node id 4
0/) ! node id 7
! Set the number of each type of element, plus the total number.
numElems=1
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/1/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_TETRA/) ! elem id 1
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(4*numElems))
elemConn=(/1,2,4,3/) ! elem id 1
else if (localPET .eq. 1) then !!! This part only for PET 1
! Set number of nodes
numNodes=4
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/2,3,5,9/)
! Allocate and fill node coordinate array.
! Since this is a 3D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/1.0,0.0,0.0, & ! node id 2
2.0,0.0,0.0, & ! node id 3
1.5,1.0,0.0, & ! node id 5
1.5,0.5,1.0/) ! node id 9
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 2
1, & ! node id 3
2, & ! node id 5
1/) ! node id 9
! Set the number of each type of element, plus the total number.
numElems=1
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/2/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_TETRA/) ! elem id 1
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(4*numElems))
elemConn=(/1,2,4,3/) ! elem id 1
else if (localPET .eq. 2) then !!! This part only for PET 2
! Set number of nodes
numNodes=4
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/2,4,5,8/)
! Allocate and fill node coordinate array.
! Since this is a 3D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/1.0,0.0,0.0, & ! node id 2
0.5,1.0,0.0, & ! node id 4
1.5,1.0,0.0, & ! node id 5
1.0,0.5,1.0/) ! node id 8
! Allocate and fill the node owner array.
! Since this Mesh is all on PET 0, it's just set to all 0.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 2
0, & ! node id 4
2, & ! node id 5
2/) ! node id 8
! Set the number of each type of element, plus the total number.
numElems=1
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/3/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_TETRA/) ! elem id 3
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(4*numElems))
elemConn=(/1,3,4,2/) ! elem id 3
else if (localPET .eq. 3) then !!! This part only for PET 3
! Set number of nodes
numNodes=4
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/4,5,6,10/)
! Allocate and fill node coordinate array.
! Since this is a 3D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(3*numNodes))
nodeCoords=(/0.5,1.0,0.0, & ! node id 4
1.5,1.0,0.0, & ! node id 5
1.0,2.0,0.0, & ! node id 6
1.0,1.5,1.0/) ! node id 10
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 4
2, & ! node id 5
3, & ! node id 6
3/) ! node id 10
! Set the number of each type of element, plus the total number.
numElems=1
! Allocate and fill the element id array.
allocate(elemIds(numElems))
elemIds=(/4/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElems))
elemTypes=(/ESMF_MESHELEMTYPE_TETRA/) ! elem id 4
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(4*numElems))
elemConn=(/1,2,4,3/) ! elem id 4
endif
endif
! Create Mesh structure in 1 step
dstMesh=ESMF_MeshCreate(parametricDim=3,spatialDim=3, &
coordSys=ESMF_COORDSYS_CART, &
nodeIds=nodeIds, nodeCoords=nodeCoords, &
nodeOwners=nodeOwners, elementIds=elemIds,&
elementTypes=elemTypes, elementConn=elemConn, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Array spec
call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=rc)
! Create dest. field
dstField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create dest. area field
dstAreaField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest_area", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create dest. frac field
dstFracField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="dest_frac", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create exact dest. field
xdstField = ESMF_FieldCreate(dstMesh, arrayspec, meshloc=ESMF_MESHLOC_ELEMENT, &
name="xdest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Init destination field to 0.0
! Should only be 1 localDE
call ESMF_FieldGet(dstField, 0, dstFarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Init destination field to 0.0
dstFarrayPtr=0.0
! Init exact destination field
! Should only be 1 localDE
call ESMF_FieldGet(xdstField, 0, xdstFarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! set interpolated function
iconn=1
do i1=1,numElems
! Loop through nodes in elem
! to compute point in center
x=0.0
y=0.0
z=0.0
do i3=1,4 ! 4 is the number of nodes in each element
inode=3*(elemConn(iconn)-1)
x=x+nodeCoords(inode+1)
y=y+nodeCoords(inode+2)
z=z+nodeCoords(inode+3)
iconn=iconn+1
enddo
x=0.25*x
y=0.25*y
z=0.25*z
! Set source function
xdstFarrayPtr(i1) = 20.0+x+y+z
enddo
! For now, Easy set interpolated function
!xdstFarrayPtr=1.0
! deallocate node data
deallocate(nodeIds)
deallocate(nodeCoords)
deallocate(nodeOwners)
! deallocate elem data
deallocate(elemIds)
deallocate(elemTypes)
deallocate(elemConn)
!!! Regrid forward from the A grid to the B grid
! Regrid store
call ESMF_FieldRegridStore( &
srcField, &
srcMaskValues=(/1/), &
dstField=dstField, &
routeHandle=routeHandle, &
regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
dstFracField=dstFracField, &
srcFracField=srcFracField, &
unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Do regrid
call ESMF_FieldRegrid(srcField, dstField, routeHandle, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldRegridRelease(routeHandle, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get the integration weights
call ESMF_FieldRegridGetArea(srcAreaField, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get the integration weights
call ESMF_FieldRegridGetArea(dstAreaField, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Check if the values are close
minerror(1) = 100000.
maxerror(1) = 0.
error = 0.
errorTot=0.0
dstmass = 0.
! get dst Field
call ESMF_FieldGet(dstField, 0, dstFarrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get exact destination Field
call ESMF_FieldGet(xdstField, 0, xdstFarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get dst area Field
call ESMF_FieldGet(dstAreaField, 0, dstAreaPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get frac Field
call ESMF_FieldGet(dstFracField, 0, dstFracptr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! 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 compute error on it.
! (Note that this is what SCRIP does)
if (dstFracptr(i1) .lt. 0.999) cycle
! write(*,*) i1,"::",dstFarrayPtr(i1),xdstFarrayPtr(i1)
if (xdstFarrayPtr(i1) .ne. 0.0) then
error=ABS(dstFarrayPtr(i1) - xdstFarrayPtr(i1))/ABS(xdstFarrayPtr(i1))
errorTot=errorTot+error
if (error > maxerror(1)) then
maxerror(1) = error
endif
if (error < minerror(1)) then
minerror(1) = error
endif
else
error=ABS(dstFarrayPtr(i1) - xdstFarrayPtr(i1))/ABS(xdstFarrayPtr(i1))
errorTot=errorTot+error
if (error > maxerror(1)) then
maxerror(1) = error
endif
if (error < minerror(1)) then
minerror(1) = error
endif
endif
enddo
srcmass(1) = 0.
! get src pointer
call ESMF_FieldGet(srcField, 0, srcFarrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get src Field
call ESMF_FieldGet(srcAreaField, 0, srcAreaptr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get frac Field
call ESMF_FieldGet(srcFracField, 0, srcFracptr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
do i1=clbnd(1),cubnd(1)
srcmass(1) = srcmass(1) + srcFracptr(i1)*srcAreaptr(i1)*srcFarrayPtr(i1)
enddo
srcmassg(1) = 0.
dstmassg(1) = 0.
call ESMF_VMAllReduce(vm, srcmass, srcmassg, 1, ESMF_REDUCE_SUM, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_VMAllReduce(vm, dstmass, dstmassg, 1, ESMF_REDUCE_SUM, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_VMAllReduce(vm, maxerror, maxerrorg, 1, ESMF_REDUCE_MAX, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_VMAllReduce(vm, minerror, minerrorg, 1, ESMF_REDUCE_MIN, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! return answer based on correct flags
csrv = .false.
if (ABS(dstmassg(1)-srcmassg(1))/srcmassg(1) < 10E-10) csrv = .true.
itrp = .false.
if (maxerrorg(1) < 10E-2) itrp = .true.
! Uncomment these calls to see some actual regrid results
if (localPet == 0) then
write(*,*) "=== 3D Cartesian Mesh ==="
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 = ", (maxerrorg(1) + minerrorg(1))/2
write(*,*) " "
endif
! Destroy the Fields
call ESMF_FieldDestroy(srcField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(dstField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(srcAreaField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(dstAreaField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(srcFracField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(dstFracField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(xdstField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Free the meshes
call ESMF_MeshDestroy(srcMesh, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_MeshDestroy(dstMesh, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! return success if we've gotten this far
rc=ESMF_SUCCESS
end subroutine test_RegridCsrv3DCartMesh