subroutine test_RegridCsrvCartPHMesh(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
integer :: iconn,inode
integer :: numQuadElems,numTriElems
integer :: numPentElems,numHexElems,numTotElems
integer :: numElemConn
! 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
! Src Mesh
!
! 2.5 8 10 --------11
! / \ / |
! 2.1 7 9 12
! | | 5 /
! | 4 | /
! | | /
! 1.0 4 ------- 5 ------- 6
! | | \ 3 |
! | 1 | \ |
! | | 2 \ |
! -0.1 1 ------- 2 ------- 3
!
! -0.1 1.0 2.1 2.5
!
! Node Id labels at corners
! Element Id labels in centers
! Setup Src Mesh
if (petCount .eq. 1) then
! Set number of nodes
numNodes=12
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/1,2,3,4,5,6,7,8,9,10,11,12/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/-0.1,-0.1, & ! node id 1
1.0,-0.1, & ! node id 2
2.1,-0.1, & ! node id 3
-0.1, 1.0, & ! node id 4
1.0, 1.0, & ! node id 5
2.1, 1.0, & ! node id 6
-0.1, 2.1, & ! node id 7
0.5, 2.5, & ! node id 8
1.0, 2.1, & ! node id 9
1.5, 2.5, & ! node id 10
2.5, 2.5, & ! node id 11
2.5, 2.1/) ! node id 12
! 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 tot and num conn.
numQuadElems=1
numTriElems=2
numPentElems=1
numHexElems=1
numTotElems=numTriElems+numQuadElems+numPentElems+numHexElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems+6*numHexElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/1,2,3,4,5/)
! Allocate and fill the element Mask array.
allocate(elemMask(numTotElems))
elemMask=(/0,1,0,0,0/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD, & ! elem id 1
ESMF_MESHELEMTYPE_TRI, & ! elem id 2
ESMF_MESHELEMTYPE_TRI, & ! elem id 3
5, & ! elem id 4
6/) ! elem id 5
! 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(numElemConn))
elemConn=(/1,2,5,4, & ! elem id 1
2,3,5, & ! elem id 2
3,6,5, & ! elem id 3
4,5,9,8,7, & ! elem id 4
5,6,12,11,10,9/) ! elem id 5
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,5/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/-0.1, -0.1, & ! node id 1
1.0, -0.1, & ! node id 2
-0.1, 1.0, & ! node id 4
1.0, 1.0 /) ! node id 5
! 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
! Set the number of each type of element, plus tot and num conn.
numQuadElems=1
numTriElems=0
numPentElems=0
numHexElems=0
numTotElems=numTriElems+numQuadElems+numPentElems+numHexElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems+6*numHexElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/1/)
! Allocate and fill the element Mask array.
allocate(elemMask(numTotElems))
elemMask=(/0/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 1
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(numElemConn))
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,6/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/1.0,-0.1, & ! node id 2
2.1,-0.1, & ! node id 3
1.0, 1.0, & ! node id 5
2.1, 1.0 /) ! node id 6
! 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
! Set the number of each type of element, plus tot and num conn.
numQuadElems=0
numTriElems=2
numPentElems=0
numHexElems=0
numTotElems=numTriElems+numQuadElems+numPentElems+numHexElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems+6*numHexElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/2,3/)
! Allocate and fill the element Mask array.
allocate(elemMask(numTotElems))
elemMask=(/1,0/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_TRI, & ! elem id 2
ESMF_MESHELEMTYPE_TRI/) ! elem id 3
! Allocate and fill the element connection type array.
allocate(elemConn(numElemConn))
elemConn=(/1,2,3, & ! elem id 2
2,4,3/) ! elem id 3
else if (localPET .eq. 2) then !!! This part only for PET 2
! Set number of nodes
numNodes=5
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/4,5,7,8,9/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/-0.1,1.0, & ! node id 4
1.0,1.0, & ! node id 5
-0.1,2.1, & ! node id 7
0.5,2.5, & ! node id 8
1.0,2.1 /) ! node id 9
! 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
2/) ! node id 9
! Set the number of each type of element, plus tot and num conn.
numQuadElems=0
numTriElems=0
numPentElems=1
numHexElems=0
numTotElems=numTriElems+numQuadElems+numPentElems+numHexElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems+6*numHexElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/4/)
! Allocate and fill the element Mask array.
allocate(elemMask(numTotElems))
elemMask=(/0/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/5/) ! elem id 4
! Allocate and fill the element connection type array.
allocate(elemConn(numElemConn))
elemConn=(/1,2,5,4,3/) ! elem id 4
else if (localPET .eq. 3) then !!! This part only for PET 3
! Set number of nodes
numNodes=6
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/5,6,9,10,11,12/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/1.0,1.0, & ! node id 5
2.1,1.0, & ! node id 6
1.0,2.1, & ! node id 9
1.5,2.5, & ! node id 10
2.5,2.5, & ! node id 11
2.5,2.1 /) ! node id 12
! Allocate and fill the node owner array.
allocate(nodeOwners(numNodes))
nodeOwners=(/0, & ! node id 5
1, & ! node id 6
2, & ! node id 9
3, & ! node id 10
3, & ! node id 11
3/) ! node id 12
! Set the number of each type of element, plus tot and num conn.
numQuadElems=0
numTriElems=0
numPentElems=0
numHexElems=1
numTotElems=numTriElems+numQuadElems+numPentElems+numHexElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems+6*numHexElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/5/)
! Allocate and fill the element Mask array.
allocate(elemMask(numTotElems))
elemMask=(/0/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/6/) ! elem id 5
! Allocate and fill the element connection type array.
allocate(elemConn(numElemConn))
elemConn=(/1,2,6,5,4,3/) ! elem id 5
endif
endif
! Create Mesh structure in 1 step
srcMesh=ESMF_MeshCreate(parametricDim=2,spatialDim=2, &
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,numTotElems
! Loop through nodes in elem
! to compute point in center
x=0.0
y=0.0
do i2=1,elemTypes(i1)
inode=2*(elemConn(iconn)-1)
x=x+nodeCoords(inode+1)
y=y+nodeCoords(inode+2)
iconn=iconn+1
enddo
x=x*(1.0/REAL(elemTypes(i1),ESMF_KIND_R8))
y=y*(1.0/REAL(elemTypes(i1),ESMF_KIND_R8))
! 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
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 !!!!!!!!!!!!!!!!!
! Dest Mesh
!
! 2.0 7 ------- 8 ------- 9
! | | \
! 1.5 | 4 | 5 10
! | | /
! 1.0 4 ------- 5 ------- 6
! | | \ 3 |
! | 1 | \ |
! | | 2 \ |
! 0.0 1 ------- 2 ------- 3
!
! 0.0 1.0 2.0 2.25
!
! Node Id labels at corners
! Element Id labels in centers
! Create Dest 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 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/0.0,0.0, & ! node id 1
1.0,0.0, & ! node id 2
2.0,0.0, & ! node id 3
0.0,1.0, & ! node id 4
1.0,1.0, & ! node id 5
2.0,1.0, & ! node id 6
0.0,2.0, & ! node id 7
1.0,2.0, & ! node id 8
2.0,2.0, & ! node id 9
2.25,1.5 /) ! 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.
numTriElems=2
numQuadElems=2
numPentElems=1
numTotElems=numTriElems+numQuadElems+numPentElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/1,2,3,4,5/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD, & ! elem id 1
ESMF_MESHELEMTYPE_TRI, & ! elem id 2
ESMF_MESHELEMTYPE_TRI, & ! elem id 3
ESMF_MESHELEMTYPE_QUAD, & ! elem id 4
5/) ! elem id 5
! 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(numElemConn))
elemConn=(/1,2,5,4, & ! elem id 1
2,3,5, & ! elem id 2
3,6,5, & ! elem id 3
4,5,8,7, & ! elem id 4
5,6,10,9,8/)! elem id 5
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,5/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/0.0,0.0, & ! node id 1
1.0,0.0, & ! node id 2
0.0,1.0, & ! node id 4
1.0,1.0 /) ! node id 5
! 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
! Set the number of each type of element, plus the total number.
numTriElems=0
numQuadElems=1
numPentElems=0
numTotElems=numTriElems+numQuadElems+numPentElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/1/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 1
! Allocate and fill the element connection type array.
! Note that entry are local indices
allocate(elemConn(numElemConn))
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,6/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/1.0,0.0, & ! node id 2
2.0,0.0, & ! node id 3
1.0,1.0, & ! node id 5
2.0,1.0 /) ! node id 6
! 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
! Set the number of each type of element, plus the total number.
numTriElems=2
numQuadElems=0
numPentElems=0
numTotElems=numTriElems+numQuadElems+numPentElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/2,3/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_TRI, & ! elem id 2
ESMF_MESHELEMTYPE_TRI/) ! elem id 3
! Allocate and fill the element connection type array.
allocate(elemConn(numElemConn))
elemConn=(/1,2,3, & ! elem id 2
2,4,3/) ! elem id 3
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=(/4,5,7,8/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/0.0,1.0, & ! node id 4
1.0,1.0, & ! node id 5
0.0,2.0, & ! node id 7
1.0,2.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 4
0, & ! node id 5
2, & ! node id 7
2/) ! node id 8
! Set the number of each type of element, plus the total number.
numTriElems=0
numQuadElems=1
numPentElems=0
numTotElems=numTriElems+numQuadElems+numPentElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/4/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numTotElems))
elemTypes=(/ESMF_MESHELEMTYPE_QUAD/) ! elem id 4
! Allocate and fill the element connection type array.
allocate(elemConn(numElemConn))
elemConn=(/1,2,4,3/) ! elem id 4
else if (localPET .eq. 3) then !!! This part only for PET 3
! Set number of nodes
numNodes=5
! Allocate and fill the node id array.
allocate(nodeIds(numNodes))
nodeIds=(/5,6,8,9,10/)
! Allocate and fill node coordinate array.
! Since this is a 2D Mesh the size is 2x the
! number of nodes.
allocate(nodeCoords(2*numNodes))
nodeCoords=(/1.0,1.0, & ! node id 5
2.0,1.0, & ! node id 6
1.0,2.0, & ! node id 8
2.0,2.0, & ! node id 9
2.25,1.5 /) ! node id 10
! 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
3/) ! node id 10
! Set the number of each type of element, plus the total number.
numTriElems=0
numQuadElems=0
numPentElems=1
numTotElems=numTriElems+numQuadElems+numPentElems
numElemConn=3*numTriElems+4*numQuadElems+ &
5*numPentElems
! Allocate and fill the element id array.
allocate(elemIds(numTotElems))
elemIds=(/5/)
! Allocate and fill the element topology type array.
allocate(elemTypes(numElemConn))
elemTypes=(/5/) ! elem id 5
! Allocate and fill the element connection type array.
allocate(elemConn(numElemConn))
elemConn=(/1,2,5,4,3/) ! elem id 5
endif
endif
! Create Mesh structure in 1 step
dstMesh=ESMF_MeshCreate(parametricDim=2,spatialDim=2, &
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,numTotElems
! Loop through nodes in elem
! to compute point in center
x=0.0
y=0.0
do i2=1,elemTypes(i1)
inode=2*(elemConn(iconn)-1)
x=x+nodeCoords(inode+1)
y=y+nodeCoords(inode+2)
iconn=iconn+1
enddo
x=x*(1.0/REAL(elemTypes(i1),ESMF_KIND_R8))
y=y*(1.0/REAL(elemTypes(i1),ESMF_KIND_R8))
! Set source function
xdstFarrayPtr(i1) = 20.0+x+y
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)
#if 0
call ESMF_MeshWrite(srcMesh,"srcMesh")
call ESMF_MeshWrite(dstMesh,"dstMesh")
#endif
!!! 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(*,*) "=== Mesh with Pentagons and Hexagons ==="
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_RegridCsrvCartPHMesh