subroutine flux_exchange(xgrid, srcField, dstField, rc)
type(ESMF_XGrid), intent(inout) :: xgrid
type(ESMF_Field), intent(inout) :: srcField(:)
type(ESMF_Field), intent(inout) :: dstField(:)
integer, intent(out), optional :: rc
integer :: localrc, i, j, nsrc, ndst, lpet, npet
type(ESMF_Field) :: f_xgrid
type(ESMF_Grid), allocatable :: srcGrid(:)
type(ESMF_Field), allocatable :: srcFrac(:), srcArea(:)
type(ESMF_Grid), allocatable :: dstGrid(:)
type(ESMF_Field), allocatable :: dstFrac(:), dstArea(:)
type(ESMF_Field), allocatable :: srcFrac2(:), 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
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
nsrc = size(srcField)
ndst = size(dstField)
allocate(srcGrid(nsrc), srcFrac(nsrc), srcFrac2(nsrc), srcArea(nsrc))
allocate(dstGrid(ndst), dstFrac(ndst), dstFrac2(ndst), dstArea(ndst))
do i = 1, size(srcField)
call ESMF_FieldGet(srcField(i), grid=srcGrid(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
srcFrac(i) = ESMF_FieldCreate(srcGrid(i), typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
srcFrac2(i) = ESMF_FieldCreate(srcGrid(i), typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
srcArea(i) = ESMF_FieldCreate(srcGrid(i), 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, size(dstField)
call ESMF_FieldGet(dstField(i), grid=dstGrid(i), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
dstFrac(i) = ESMF_FieldCreate(dstGrid(i), typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
dstFrac2(i) = ESMF_FieldCreate(dstGrid(i), 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), 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
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
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
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
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((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
enddo
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)
call ESMF_XGridDestroy(xgrid,rc=localrc)
if(present(rc)) rc = ESMF_SUCCESS
end subroutine flux_exchange