subroutine test_regridMeshToGridSph3D(rc)
integer, intent(out) :: rc
logical :: correct
integer :: localrc
type(ESMF_Mesh) :: srcMesh
type(ESMF_Grid) :: dstGrid
type(ESMF_Field) :: srcField
type(ESMF_Field) :: dstField
type(ESMF_Field) :: xdstField
type(ESMF_Array) :: dstArray
type(ESMF_Array) :: xdstArray
type(ESMF_Array) :: lonArrayA
type(ESMF_Array) :: srcArrayA
type(ESMF_RouteHandle) :: routeHandle
type(ESMF_ArraySpec) :: arrayspec
type(ESMF_VM) :: vm
real(ESMF_KIND_R8), pointer :: farrayPtrXC(:,:,:), farrayPtr1D(:)
real(ESMF_KIND_R8), pointer :: farrayPtrYC(:,:,:),farrayPtrZC(:,:,:)
real(ESMF_KIND_R8), pointer :: farrayPtr(:,:,:),farrayPtr2(:,:,:)
real(ESMF_KIND_R8), pointer :: xfarrayPtr(:,:,:)
integer :: clbnd(3),cubnd(3)
integer :: fclbnd(3),fcubnd(3)
integer :: i1,i2,i3, index(3)
integer :: lDE, localDECount
real(ESMF_KIND_R8) :: coord(3)
character(len=ESMF_MAXSTR) :: string
integer dst_nx,dst_ny,dst_nz
integer num_arrays
real(ESMF_KIND_R8) :: dx,dy,dz
real(ESMF_KIND_R8) :: dst_minx,dst_miny,dst_minz
real(ESMF_KIND_R8) :: dst_maxx,dst_maxy,dst_maxz
real(ESMF_KIND_R8) :: x,y,z
real(ESMF_KIND_R8) :: lat_rad, lon_rad
real(ESMF_KIND_R8) :: err, relErr, maxRelErr
integer :: spherical_grid
integer, pointer :: larrayList(:)
integer :: localPet, petCount
integer, pointer :: nodeIds(:),nodeOwners(:)
real(ESMF_KIND_R8), pointer :: nodeCoords(:)
integer, pointer :: elemIds(:),elemTypes(:),elemConn(:)
integer :: numNodes, numElems
! result code
integer :: finalrc
! degree to rad conversion
real(ESMF_KIND_R8),parameter :: &
DEG2RAD = 3.141592653589793_ESMF_KIND_R8/180.0_ESMF_KIND_R8
! 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
! If we don't have 1 or 4 PETS then exit successfully
if ((petCount .ne. 1) .and. (petCount .ne. 4)) then
rc=ESMF_SUCCESS
return
endif
! Establish the resolution of the grids
dst_nx = 10
dst_ny = 10
dst_nz = 5
! Establish the coordinates of the grids
dst_minx = 1.1
dst_miny = 1.1
dst_minz = 1.1
dst_maxx = 19.9
dst_maxy = 19.9
dst_maxz = 1.9
! 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=(/1.0,1.0,1.0, & ! node id 1
10.0,1.0,1.0, & ! node id 2
20.0,1.0,1.0, & ! node id 3
1.0,10.0,1.0, & ! node id 4
10.0,10.0,1.0, & ! node id 5
20.0,10.0,1.0, & ! node id 6
1.0,20.0,1.0, & ! node id 7
10.0,20.0,1.0, & ! node id 8
20.0,20.0,1.0, & ! node id 9
1.0,1.0,2.0, & ! node id 10
10.0,1.0,2.0, & ! node id 11
20.0,1.0,2.0, & ! node id 12
1.0,10.0,2.0, & ! node id 13
10.0,10.0,2.0, & ! node id 14
20.0,10.0,2.0, & ! node id 15
1.0,20.0,2.0, & ! node id 16
10.0,20.0,2.0, & ! node id 17
20.0,20.0,2.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 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=(/1.0,1.0,1.0, & ! node id 1
10.0,1.0,1.0, & ! node id 2
1.0,10.0,1.0, & ! node id 4
10.0,10.0,1.0, & ! node id 5
1.0,1.0,2.0, & ! node id 10
10.0,1.0,2.0, & ! node id 11
1.0,10.0,2.0, & ! node id 13
10.0,10.0,2.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 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=(/10.0,1.0,1.0, & ! node id 2
20.0,1.0,1.0, & ! node id 3
10.0,10.0,1.0, & ! node id 5
20.0,10.0,1.0, & ! node id 6
10.0,1.0,2.0, & ! node id 11
20.0,1.0,2.0, & ! node id 12
10.0,10.0,2.0, & ! node id 14
20.0,10.0,2.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 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=(/1.0,10.0,1.0, & ! node id 4
10.0,10.0,1.0, & ! node id 5
1.0,20.0,1.0, & ! node id 7
10.0,20.0,1.0, & ! node id 8
1.0,10.0,2.0, & ! node id 13
10.0,10.0,2.0, & ! node id 14
1.0,20.0,2.0, & ! node id 16
10.0,20.0,2.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 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=(/10.0,10.0,1.0, & ! node id 5
20.0,10.0,1.0, & ! node id 6
10.0,20.0,1.0, & ! node id 8
20.0,20.0,1.0, & ! node id 9
10.0,10.0,2.0, & ! node id 14
20.0,10.0,2.0, & ! node id 15
10.0,20.0,2.0, & ! node id 17
20.0,20.0,2.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 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, &
nodeIds=nodeIds, nodeCoords=nodeCoords, &
nodeOwners=nodeOwners, elementIds=elemIds,&
elementTypes=elemTypes, elementConn=elemConn, &
coordSys=ESMF_COORDSYS_SPH_DEG, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source field
call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=rc)
srcField = ESMF_FieldCreate(srcMesh, arrayspec, &
name="source", 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, farrayPtr1D, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! set interpolated function
i2=1
do i1=1,numNodes
if (nodeOwners(i1) .eq. localPet) then
! Get coordinates
x=nodeCoords(3*i1-2)
y=nodeCoords(3*i1-1)
z=nodeCoords(3*i1)
! Compute Lat/Lon in rad
lon_rad = DEG2RAD*(x)
lat_rad = DEG2RAD*(90.-y)
! Set source function
!farrayPtr1D(i2) = 20.0+x+y+z
!farrayPtr1D(i2) = z*(2. + cos(lon_rad)**2.*cos(2.*lat_rad))
farrayPtr1D(i2) = z*5.0+2.0+ cos(lon_rad)**2.*cos(2.*lat_rad)
!farrayPtr1D(i2) = 1.0+z*0.01+ cos(lat_rad)**2.*cos(2.*lon_rad)
!farrayPtr1D(i2) = z
! Advance to next owner
i2=i2+1
endif
enddo
! deallocate node data
deallocate(nodeIds)
deallocate(nodeCoords)
deallocate(nodeOwners)
! deallocate elem data
deallocate(elemIds)
deallocate(elemTypes)
deallocate(elemConn)
! setup dest. grid
dstGrid=ESMF_GridCreateNoPeriDim(minIndex=(/1,1,1/),maxIndex=(/dst_nx,dst_ny,dst_nz/), &
coordSys=ESMF_COORDSYS_SPH_DEG, regDecomp=(/2,2,1/), indexflag=ESMF_INDEX_GLOBAL, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source/destination fields
call ESMF_ArraySpecSet(arrayspec, 3, ESMF_TYPEKIND_R8, rc=rc)
dstField = ESMF_FieldCreate(dstGrid, arrayspec, &
staggerloc=ESMF_STAGGERLOC_CENTER_VCENTER, name="dest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Exact destination field
xdstField = ESMF_FieldCreate(dstGrid, arrayspec, &
staggerloc=ESMF_STAGGERLOC_CENTER_VCENTER, name="dest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridAddCoord(dstGrid, staggerloc=ESMF_STAGGERLOC_CENTER_VCENTER, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get number of local DEs
call ESMF_GridGet(dstGrid, localDECount=localDECount, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get arrays
! dstArray
call ESMF_FieldGet(dstField, array=dstArray, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(xdstField, array=xdstArray, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! srcArrayA
call ESMF_FieldGet(srcField, array=srcArrayA, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Destination grid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Get memory and set coords for dst
do lDE=0,localDECount-1
!! get coords
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=1, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrXC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=2, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrYC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=3, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrZC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(dstField, lDE, farrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(xdstField, lDE, xfarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!! set coords
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
do i3=clbnd(3),cubnd(3)
! Set source coordinates
farrayPtrXC(i1,i2,i3) = ((dst_maxx-dst_minx)*REAL(i1-1)/REAL(dst_nx-1))+dst_minx
farrayPtrYC(i1,i2,i3) = ((dst_maxy-dst_miny)*REAL(i2-1)/REAL(dst_ny-1))+dst_miny
farrayPtrZC(i1,i2,i3) = ((dst_maxz-dst_minz)*REAL(i3-1)/REAL(dst_nz-1))+dst_minz
! Put in more convenient form
x=farrayPtrXC(i1,i2,i3)
y=farrayPtrYC(i1,i2,i3)
z=farrayPtrZC(i1,i2,i3)
! Compute Lat/Lon in rad
lon_rad = DEG2RAD*(x)
lat_rad = DEG2RAD*(90.-y)
! initialize exact destination field
!xfarrayPtr(i1,i2,i3)=z*(2. + cos(lon_rad)**2.*cos(2.*lat_rad))
xfarrayPtr(i1,i2,i3)= z*5.0+2.0+ cos(lon_rad)**2.*cos(2.*lat_rad)
!xfarrayPtr(i1,i2,i3)= 1.0+z*0.01+ cos(lat_rad)**2.*cos(2.*lon_rad)
!xfarrayPtr(i1,i2,i3)=z
! Initialize destination field
farrayPtr(i1,i2,i3)=0.0
enddo
enddo
enddo
enddo ! lDE
!!! Regrid forward from the A grid to the B grid
! Regrid store
call ESMF_FieldRegridStore( &
srcField, &
dstField=dstField, &
routeHandle=routeHandle, &
regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
lineType=ESMF_LINETYPE_GREAT_CIRCLE, &
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
! Init max
maxRelErr=0.0
! Check error
do lDE=0,localDECount-1
!! get coords
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, &
farrayPtr=farrayPtrXC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, &
farrayPtr=farrayPtrYC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=3, &
farrayPtr=farrayPtrZC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(dstField, lDE, farrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(xdstField, lDE, xfarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!! check error
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
do i3=clbnd(3),cubnd(3)
! Calculate error
err=abs(farrayPtr(i1,i2,i3)-xfarrayPtr(i1,i2,i3))
! Calculate relative error
if (xfarrayPtr(i1,i2,i3) .ne. 0.0) then
relErr=err/xfarrayPtr(i1,i2,i3)
else
relErr=err
endif
!! Calculate max
if (relErr > maxRelErr) maxRelErr=relErr
!! if error is too big report an error
if (relErr > 0.009) then
correct=.false.
endif
enddo
enddo
enddo
enddo ! lDE
! write(*,*) "Max Rel. Error= ",maxRelErr
#if 0
call ESMF_MeshWrite(srcMesh,"srcMesh",rc=localrc)
call ESMF_GridWriteVTK(dstGrid,staggerloc=ESMF_STAGGERLOC_CENTER, &
filename="dstGrid", array1=dstArray, array2=xdstArray, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
#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
! Free the grids
call ESMF_MeshDestroy(srcMesh, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridDestroy(dstGrid, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! return answer based on correct flag
if (correct) then
rc=ESMF_SUCCESS
else
rc=ESMF_FAILURE
endif
end subroutine test_regridMeshToGridSph3D