subroutine flux_exchange_sph_mesh(xgrid, scheme, rc)
type(ESMF_XGrid), intent(inout) :: xgrid
integer, intent(in), optional :: scheme
integer, intent(out), optional :: rc
integer :: localrc, i, j, nsrc, ndst, lpet, npet
type(ESMF_Field) :: f_xgrid
type(ESMF_Mesh), allocatable :: srcGrid(:)
type(ESMF_Field), allocatable :: srcFrac(:), srcArea(:)
type(ESMF_Mesh), allocatable :: dstGrid(:)
type(ESMF_Field), allocatable :: dstFrac(:), dstArea(:)
type(ESMF_Field), allocatable :: srcFrac2(:)
type(ESMF_Field), allocatable :: dstFrac2(:)
type(ESMF_RouteHandle), allocatable :: s2d_rh(:,:)
type(ESMF_RouteHandle), allocatable :: d2s_rh(:,:)
type(ESMF_RouteHandle), allocatable :: s2x_rh(:), x2s_rh(:)
type(ESMF_RouteHandle), allocatable :: d2x_rh(:), x2d_rh(:)
real(ESMF_KIND_R8), pointer :: src(:), dst(:), exf(:)
real(ESMF_KIND_R8), pointer :: src_area(:), dst_area(:), exf_area(:)
real(ESMF_KIND_R8), pointer :: src_frac(:), dst_frac(:), exf_frac(:)
real(ESMF_KIND_R8), pointer :: src_frac2(:), dst_frac2(:)
real(ESMF_KIND_R8) :: srcsum(3), allsrcsum(3), scale=2.0, exf_tarea, exf_tflux
type(ESMF_VM) :: vm
type(ESMF_Field), allocatable :: srcField(:)
type(ESMF_Field), allocatable :: dstField(:)
integer :: l_scheme
integer :: sideAGC, sideBGC, sideAMC, sideBMC
real(ESMF_KIND_R8) :: global_sum
character(len=1) :: cids(10) = (/'0','1','2','3','4','5','6','7','8','9'/)
l_scheme = ESMF_REGRID_SCHEME_REGION3D
if(present(scheme)) l_scheme = scheme
call ESMF_VMGetCurrent(vm=vm, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMGet(vm, localpet=lpet, petCount=npet, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!------------------------------------
! build Fields on the Grids
!------------------------------------
! create a Field on the xgrid
f_xgrid = ESMF_FieldCreate(xgrid=xgrid, TYPEKIND=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldGet(f_xgrid, farrayPtr=exf, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_XGridGet(xgrid, &
sideAGridCount=sideAGC, sideBGridCount=sideBGC, &
sideAMeshCount=sideAMC, sideBMeshCount=sideBMC, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
nsrc = sideAMC
ndst = sideBMC
allocate(srcGrid(nsrc), srcField(nsrc), srcFrac(nsrc), srcFrac2(nsrc), srcArea(nsrc))
allocate(dstGrid(ndst), dstField(ndst), dstFrac(ndst), dstFrac2(ndst), dstArea(ndst))
call ESMF_XGridGet(xgrid, sideAMesh=srcGrid, sideBMesh=dstGrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i = 1, nsrc
srcField(i) = ESMF_FieldCreate(srcGrid(i), meshloc=ESMF_MESHLOC_ELEMENT, &
typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
srcFrac(i) = ESMF_FieldCreate(srcGrid(i), meshloc=ESMF_MESHLOC_ELEMENT, &
typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
srcFrac2(i) = ESMF_FieldCreate(srcGrid(i), meshloc=ESMF_MESHLOC_ELEMENT, &
typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
srcArea(i) = ESMF_FieldCreate(srcGrid(i), meshloc=ESMF_MESHLOC_ELEMENT, &
typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegridGetArea(srcArea(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
enddo
do i = 1, ndst
dstField(i) = ESMF_FieldCreate(dstGrid(i), meshloc=ESMF_MESHLOC_ELEMENT, &
typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
dstFrac(i) = ESMF_FieldCreate(dstGrid(i), meshloc=ESMF_MESHLOC_ELEMENT, &
typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
dstFrac2(i) = ESMF_FieldCreate(dstGrid(i), meshloc=ESMF_MESHLOC_ELEMENT, &
typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldGet(dstFrac(i), localDe=0, farrayPtr=dst_frac, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
dstArea(i) = ESMF_FieldCreate(dstGrid(i), meshloc=ESMF_MESHLOC_ELEMENT, &
typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegridGetArea(dstArea(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
enddo
allocate(s2d_rh(size(srcField), size(dstField)), d2s_rh(size(dstField), size(srcField)))
allocate(s2x_rh(size(srcField)), x2s_rh(size(srcField)))
allocate(d2x_rh(size(dstField)), x2d_rh(size(dstField)))
do i = 1, size(srcField)
do j = 1, size(dstField)
call ESMF_FieldRegridStore(srcField=srcField(i), dstField=dstField(j), &
regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
routehandle=s2d_rh(i,j), &
unmappedaction = ESMF_UNMAPPEDACTION_IGNORE, &
srcFracField=srcFrac(i), dstFracField=dstFrac(j), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegridStore(srcField=dstField(j), dstField=srcField(i), &
regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
routehandle=d2s_rh(j,i), &
unmappedaction = ESMF_UNMAPPEDACTION_IGNORE, &
srcFracField=dstFrac(j), dstFracField=srcFrac(i), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
enddo
enddo
do i = 1, size(srcField)
call ESMF_FieldRegridStore(xgrid, srcField=srcField(i), dstField=f_xgrid, &
routehandle=s2x_rh(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegridStore(xgrid, srcField=f_xgrid, dstField=srcField(i), &
routehandle=x2s_rh(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!call ESMF_MeshWrite(srcGrid(i), cids(i)//'_srcmesh.vtk', rc=localrc)
!if (ESMF_LogFoundError(localrc, &
! ESMF_ERR_PASSTHRU, &
! ESMF_CONTEXT, rcToReturn=rc)) return
enddo
do i = 1, size(dstField)
call ESMF_FieldRegridStore(xgrid, srcField=dstField(i), dstField=f_xgrid, &
routehandle=d2x_rh(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldRegridStore(xgrid, srcField=f_xgrid, dstField=dstField(i), &
routehandle=x2d_rh(i), dstFracField=dstFrac(i), dstMergeFracField=dstFrac2(i), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!call ESMF_MeshWrite(dstGrid(i), cids(i)//'_dstmesh.vtk', rc=localrc)
!if (ESMF_LogFoundError(localrc, &
! ESMF_ERR_PASSTHRU, &
! ESMF_CONTEXT, rcToReturn=rc)) return
enddo
!----------------------------------------------------
! Compute flux integrals
! Initialize src flux to constant
!----------------------------------------------------
exf = 0.
do i = 1, size(srcField)
call ESMF_FieldGet(srcField(i), farrayPtr=src, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
src = scale
enddo
! Perform flux exchange
do i = 1, size(srcField)
call ESMF_FieldRegrid(srcField=srcField(i), dstField=f_xgrid, &
routehandle=s2x_rh(i), zeroregion=ESMF_REGION_EMPTY, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
enddo
! make sure flux is conserved on XGrid
!call ESMF_MeshWrite(xgrid%xgtypep%mesh, 'xgrid.vtk', rc=localrc)
!if (ESMF_LogFoundError(localrc, &
! ESMF_ERR_PASSTHRU, &
! ESMF_CONTEXT, rcToReturn=rc)) return
allocate(exf_area(lbound(exf,1):ubound(exf,1)))
allocate(exf_frac(lbound(exf,1):ubound(exf,1)))
call ESMF_XGridGet(xgrid, area=exf_area, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
exf_frac = 1.0
call compute_flux1D(vm, exf, exf_area, exf_frac, allsrcsum, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
if(lpet == 0) print *, ' xgrid flux and area: ', allsrcsum
!if(abs(allsrcsum(1) - allsrcsum(2)*scale) .gt. 1.e-10) then
! call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, &
! msg="- inconsistent flux and area found", &
! ESMF_CONTEXT, rcToReturn=rc)
! return
!endif
exf_tflux = allsrcsum(1)
exf_tarea = allsrcsum(2)
deallocate(exf_area, exf_frac)
!make sure flux is conserved on dst Fields
global_sum = 0.
do i = 1, size(dstField)
call ESMF_FieldRegrid(srcField=f_xgrid, dstField=dstField(i), &
routehandle=x2d_rh(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldGet(dstField(i), farrayPtr=dst, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! fraction
call ESMF_FieldGet(dstFrac(i), farrayPtr=dst_frac, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldGet(dstFrac2(i), farrayPtr=dst_frac2, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! area
call ESMF_FieldGet(dstArea(i), farrayPtr=dst_area, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!call compute_flux2D(vm, dst, dst_area, dst_frac, dst_frac2, allsrcsum, dstflux=.true., rc=localrc)
!if (ESMF_LogFoundError(localrc, &
! ESMF_ERR_PASSTHRU, &
! ESMF_CONTEXT, rcToReturn=rc)) return
!if(lpet == 0) print *, 'dst flux and area: ', allsrcsum
!if(ndst == 1) then
! if((abs(exf_tarea - allsrcsum(2)) .gt. 1.e-10) .or. &
! (abs(exf_tflux - allsrcsum(1)) .gt. 1.e-10)) then
! call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, &
! msg="- inconsistent flux and area found", &
! ESMF_CONTEXT, rcToReturn=rc)
! return
! endif
!else
! call compute_flux2D(vm, dst, dst_area, dst_frac, dst_frac2, allsrcsum, dstflux=.true., rc=localrc)
! if (ESMF_LogFoundError(localrc, &
! ESMF_ERR_PASSTHRU, &
! ESMF_CONTEXT, rcToReturn=rc)) return
! if(lpet == 0) print *, 'dst flux and area using frac2: ', allsrcsum
! global_sum = global_sum + allsrcsum(1)
!endif
enddo
! make sure going to multiple Grids also conserve global flux
!if(ndst .gt. 1) then
! if ((abs(exf_tflux - global_sum) .gt. 1.e-10)) then
! call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_WRONG, &
! msg="- inconsistent flux and area found", &
! ESMF_CONTEXT, rcToReturn=rc)
! return
! endif
!endif
do i = 1, size(dstField)
call ESMF_FieldRegrid(srcField=dstField(i), dstField=f_xgrid, &
routehandle=d2x_rh(i), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
enddo
do i = 1, size(srcField)
call ESMF_FieldRegrid(srcField=f_xgrid, dstField=srcField(i), &
routehandle=x2s_rh(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_FieldGet(srcField(i), farrayPtr=src, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
enddo
!----------------------------------------------------
! clean up
!----------------------------------------------------
do i = 1, size(srcField)
call ESMF_FieldDestroy(srcField(i), rc=localrc)
call ESMF_FieldDestroy(srcArea(i), rc=localrc)
call ESMF_FieldDestroy(srcFrac(i), rc=localrc)
call ESMF_FieldDestroy(srcFrac2(i), rc=localrc)
call ESMF_RoutehandleRelease(s2x_rh(i), rc=localrc)
call ESMF_RoutehandleRelease(x2s_rh(i), rc=localrc)
enddo
do i = 1, size(dstField)
call ESMF_FieldDestroy(dstField(i), rc=localrc)
call ESMF_FieldDestroy(dstArea(i), rc=localrc)
call ESMF_FieldDestroy(dstFrac(i), rc=localrc)
call ESMF_FieldDestroy(dstFrac2(i), rc=localrc)
call ESMF_RoutehandleRelease(d2x_rh(i), rc=localrc)
call ESMF_RoutehandleRelease(x2d_rh(i), rc=localrc)
enddo
do i = 1, size(srcField)
do j = 1, size(dstField)
call ESMF_RoutehandleRelease(s2d_rh(i,j), rc=localrc)
call ESMF_RoutehandleRelease(d2s_rh(j,i), rc=localrc)
enddo
enddo
deallocate(srcArea, srcFrac, dstArea, dstFrac)
deallocate(s2d_rh, d2s_rh)
deallocate(s2x_rh, x2s_rh)
deallocate(d2x_rh, x2d_rh)
! Destroy xgrid
call ESMF_XGridDestroy(xgrid,rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! If rc present, return success
if(present(rc)) rc = ESMF_SUCCESS
end subroutine flux_exchange_sph_mesh