program ESMF_XGridMaskingUTest
!------------------------------------------------------------------------------
#include "ESMF.h"
!==============================================================================
!BOPI
! !PROGRAM: ESMF_XGridMaskingUTest - Unit tests for Field Create and Get methods
!
! !DESCRIPTION:
!
! The code in this file drives F90 Field Create and Get unit tests.
! The companion folder Fieldsrc contains the definitions for the
! Field methods.
!EOPI
!-----------------------------------------------------------------------------
! !USES:
use ESMF_TestMod ! test methods
use ESMF
implicit none
!------------------------------------------------------------------------------
! The following line turns the CVS identifier string into a printable variable.
character(*), parameter :: version = &
'$Id$'
! cumulative result: count failures; no failures equals "all pass"
integer :: result = 0
! individual test result code
integer :: rc = 1
! individual test failure message
character(ESMF_MAXSTR) :: failMsg
character(512) :: name
call ESMF_TestStart(ESMF_SRCLINE, rc=rc)
if (rc /= ESMF_SUCCESS) call ESMF_Finalize(endflag=ESMF_END_ABORT)
!------------------------------------------------------------------------
!NEX_UTest
! Create an XGrid in 2D
print *, 'Starting test2d_masking_single'
call test2d_masking_single(rc)
write(failMsg, *) ""
write(name, *) "Creating an XGrid from 2 adjacent Grids with masking"
call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)
!------------------------------------------------------------------------
!NEX_UTest
! Create an XGrid in 3D
print *, 'Starting test3d_masking_regional'
call test3d_masking_regional(rc)
write(failMsg, *) ""
write(name, *) "Creating an XGrid in 3D with Grid masking"
call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)
!------------------------------------------------------------------------
!NEX_UTest
! Create an XGrid in 3D
print *, 'Starting test3d_masking_global'
call test3d_masking_global(rc)
write(failMsg, *) ""
write(name, *) "Creating an XGrid in 3D with global Grid masking"
call ESMF_Test((rc.eq.ESMF_SUCCESS), name, failMsg, result, ESMF_SRCLINE)
call ESMF_TestEnd(ESMF_SRCLINE)
contains
#define ESMF_METHOD "ESMF_TESTS"
!------------------------------------------------------------------------
subroutine test2d_masking_single(rc)
integer, intent(out) :: rc
integer :: localrc, npet, lpet
type(ESMF_XGrid) :: xgrid
type(ESMF_Grid) :: sideA(1), sideB(1)
type(ESMF_VM) :: vm
type(ESMF_Field) :: srcField(2), dstField(1)
rc = ESMF_SUCCESS
localrc = ESMF_SUCCESS
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
! partially overlap
sideA(1) = make_grid(4,4,1.,1.,0.,0.,&
msx=3., mex=5., msy=0., mey=5., maskvalue=2, field=srcField(1),rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
sideB(1) = make_grid(4,8,1.,0.5,0.,0.,&
msx=0., mex=1., msy=0., mey=5., maskvalue=2, field=dstField(1),rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!xgrid = ESMF_XGridCreate(sideA, sideB, rc=localrc)
!if (ESMF_LogFoundError(localrc, &
! ESMF_ERR_PASSTHRU, &
! ESMF_CONTEXT, rcToReturn=rc)) return
xgrid = ESMF_XGridCreate(sideAGrid=sideA, sideBGrid=sideB, &
sideAMaskValues=(/2/), sideBMaskValues=(/2/), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange(xgrid, srcField(1:1), dstField(1:1), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
xgrid = ESMF_XGridCreate( &
sideAGrid=(/make_grid(4,4,1.,1.,0.,0.,msx=3.,mex=4., msy=0., mey=4., &
maskvalue=2, field=srcField(1), rc=localrc), &
make_grid(4,4,1.,1.,0.,0.,msx=0.,mex=1., msy=0., mey=4., &
maskvalue=2, field=srcField(2), rc=localrc)/), &
sideBGrid=(/make_grid(8,8,1.,1.,0.,0.,field=dstField(1), rc=localrc)/), &
sideAMaskValues=(/2,3,4/), sideBMaskValues=(/2,3,4/), &
rc=localrc)
call flux_exchange(xgrid, srcField(1:2), dstField(1:1), rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
end subroutine test2d_masking_single
subroutine test3d_masking_regional(rc)
integer, intent(out) :: rc
integer :: localrc, npet
type(ESMF_XGrid) :: xgrid
type(ESMF_Grid) :: sideA(2), sideB(1)
type(ESMF_VM) :: vm
rc = ESMF_SUCCESS
localrc = ESMF_SUCCESS
call ESMF_VMGetCurrent(vm=vm, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMGet(vm, petcount=npet, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! partially overlap
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(4,4,1.,1.,0.,0.,msx=3.,mex=4., msy=0., mey=4., maskvalue=2, rc=localrc), &
make_grid_sph(4,4,1.,1.,0.,0.,msx=0.,mex=1., msy=0., mey=4., maskvalue=2, rc=localrc)/), &
sideBGrid=(/make_grid_sph(8,8,1.,1.,0.,0.,rc=localrc)/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange_sph(xgrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Grids near north pole
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(4,4,1.,1.,0.,86.,msx=3.,mex=4., msy=-90., mey=90., maskvalue=2, rc=localrc), &
make_grid_sph(4,4,1.,1.,0.,86.,msx=0.,mex=1., msy=-90., mey=90., maskvalue=2, rc=localrc)/), &
sideBGrid=(/make_grid_sph(8,8,1.,1.,0.,82.,rc=localrc)/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange_sph(xgrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Grids near south pole
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(4,4,1.,1.,0.,-90.,msx=3.,mex=4., msy=-90., mey=90., maskvalue=2, rc=localrc), &
make_grid_sph(4,4,1.,1.,0.,-90.,msx=0.,mex=1., msy=-90., mey=90., maskvalue=2, rc=localrc)/), &
sideBGrid=(/make_grid_sph(8,8,1.,1.,0.,-90.,rc=localrc)/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange_sph(xgrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(40,40,1.,1.,0.,0.,msx=30.,mex=40., msy=0., mey=40., maskvalue=2, rc=localrc), &
make_grid_sph(40,40,1.,1.,0.,0.,msx=0.,mex=10., msy=0., mey=40., maskvalue=2, rc=localrc)/), &
sideBGrid=(/make_grid_sph(80,80,1.,1.,0.,0.,rc=localrc)/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(40,40,1.,1.,0.,0.,rc=localrc), &
make_grid_sph(40,40,1.,1.,20.5,30.5,rc=localrc)/), &
sideBGrid=(/make_grid_sph(80,80,1.,1.,0.,0.,rc=localrc)/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! With masking in XGrid
! partially overlap
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(4,4,1.,1.,0.,0.,msx=3.,mex=4., msy=0., mey=4., maskvalue=2, rc=localrc), &
make_grid_sph(4,4,1.,1.,0.,0.,msx=0.,mex=1., msy=0., mey=4., maskvalue=2, rc=localrc)/), &
sideBGrid=(/make_grid_sph(8,8,1.,1.,0.,0.,rc=localrc)/), &
sideAMaskValues=(/2,3,4/), sideBMaskValues=(/2,3,4/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange_sph(xgrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Grids near south pole
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(4,4,1.,1.,0.,-90.,msx=3.,mex=4., msy=-90., mey=90., maskvalue=2, rc=localrc), &
make_grid_sph(4,4,1.,1.,0.,-90.,msx=0.,mex=1., msy=-90., mey=90., maskvalue=2, rc=localrc)/), &
sideBGrid=(/make_grid_sph(8,8,1.,1.,0.,-90.,rc=localrc)/), &
sideAMaskValues=(/2,3,4/), sideBMaskValues=(/2,3,4/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange_sph(xgrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(40,40,1.,1.,0.,0.,msx=30.,mex=40., msy=0., mey=40., maskvalue=2, rc=localrc), &
make_grid_sph(40,40,1.,1.,0.,0.,msx=0.,mex=10., msy=0., mey=40., maskvalue=2, rc=localrc)/), &
sideBGrid=(/make_grid_sph(80,80,1.,1.,0.,0.,rc=localrc)/), &
sideAMaskValues=(/2,3,4/), sideBMaskValues=(/2,3,4/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
xgrid = ESMF_XGridCreate(&
sideAGrid=(/make_grid_sph(40,40,1.,1.,0.,0.,rc=localrc), &
make_grid_sph(40,40,1.,1.,20.5,30.5,rc=localrc)/), &
sideBGrid=(/make_grid_sph(80,80,1.,1.,0.,0.,rc=localrc)/), &
sideAMaskValues=(/2,3,4/), sideBMaskValues=(/2,3,4/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!xgrid = ESMF_XGridCreate((/make_grid_sph(40,90,1.,2.,0.,-90.,rc=localrc), &
! make_grid_sph(40,60,1.,3.,20.5,-90.,rc=localrc)/), &
! (/make_grid_sph(80,80,1.,1.,0.,0.,rc=localrc)/), &
! rc=localrc)
!if (ESMF_LogFoundError(localrc, &
! ESMF_ERR_PASSTHRU, &
! ESMF_CONTEXT, rcToReturn=rc)) return
!call flux_exchange_sph(xgrid, rc=localrc)
!if (ESMF_LogFoundError(localrc, &
! ESMF_ERR_PASSTHRU, &
! ESMF_CONTEXT, rcToReturn=rc)) return
end subroutine test3d_masking_regional
subroutine test3d_masking_global(rc)
integer, intent(out) :: rc
integer :: localrc, npet
type(ESMF_XGrid) :: xgrid
type(ESMF_Grid) :: sideA(2), sideB(1)
type(ESMF_VM) :: vm
rc = ESMF_SUCCESS
localrc = ESMF_SUCCESS
call ESMF_VMGetCurrent(vm=vm, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMGet(vm, petcount=npet, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
xgrid = ESMF_XGridCreate(&
sideAGrid=(/ &
make_grid_sph(12,9,30.,20.,0.,-90.,msx=0., mex=160., msy=-90., mey=90., &
maskvalue=3, scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc), &
make_grid_sph(9,6,40.,30.,0.,-90., msx=200., mex=360., msy=-90., mey=90., &
maskvalue=4, scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc)/), &
sideBGrid=(/make_grid_sph(12,18,30.,10.,0.,-90.,scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc)/), &
sideAMaskValues=(/2,3,4/), sideBMaskValues=(/2,3,4/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange_sph(xgrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
xgrid = ESMF_XGridCreate(&
sideAGrid=(/ &
make_grid_sph(12,9,30.,20.,0.,-90.,msx=0., mex=160., msy=-90., mey=90., &
maskvalue=3, scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc), &
make_grid_sph(9,9,40.,20.,0.,-90., msx=200., mex=360., msy=-90., mey=90., &
maskvalue=4, scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc)/), &
sideBGrid=(/make_grid_sph(12,18,30.,10.,0.,-90.,scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc)/), &
sideAMaskValues=(/2,3,4/), sideBMaskValues=(/2,3,4/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange_sph(xgrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! partially overlap
xgrid = ESMF_XGridCreate(&
sideAGrid=(/ &
make_grid_sph(120,90,3.,2.,0.,-90.,msx=0., mex=160., msy=-90., mey=90., &
maskvalue=3, scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc), &
make_grid_sph(45,120,8.,1.5,0.,-90., msx=200., mex=360., msy=-90., mey=90., &
maskvalue=4, scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc)/), &
sideBGrid=(/make_grid_sph(120,120,3.,1.5,0.,-90.,scheme=ESMF_REGRID_SCHEME_FULL3D,rc=localrc)/), &
sideAMaskValues=(/2,3,4/), sideBMaskValues=(/2,3,4/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call flux_exchange_sph(xgrid, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
end subroutine test3d_masking_global
!------------------------------------------------------------------------
! Utility methods
!------------------------------------------------------------------------
function make_grid(atm_nx, atm_ny, atm_dx, atm_dy, atm_sx, atm_sy, &
msx, msy, mex, mey, maskvalue, field, rc)
! return value
type(ESMF_Grid) :: make_grid
! arguments
integer, intent(in) :: atm_nx, atm_ny
real(ESMF_KIND_R4), intent(in) :: atm_dx, atm_dy
real(ESMF_KIND_R4), intent(in) :: atm_sx, atm_sy
real(ESMF_KIND_R4), intent(in), optional :: msx, msy, mex, mey
integer(ESMF_KIND_I4), intent(in), optional :: maskvalue
type(ESMF_Field), intent(out), optional :: field
integer, intent(out), optional :: rc
! local variables
integer :: localrc, i, j
real(ESMF_KIND_R8), pointer :: coordX(:), coordY(:)
integer(ESMF_KIND_I4), pointer :: maskptr(:,:)
real(ESMF_KIND_R8) :: startx, starty
make_grid = ESMF_GridCreateNoPeriDim(maxIndex=(/atm_nx, atm_ny/), &
coordSys=ESMF_COORDSYS_CART, &
indexflag=ESMF_INDEX_GLOBAL, &
coordDep1=(/1/), &
coordDep2=(/2/), &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridAddItem(make_grid, staggerloc=ESMF_STAGGERLOC_CENTER, &
itemflag=ESMF_GRIDITEM_MASK, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetItem(make_grid, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
itemflag=ESMF_GRIDITEM_MASK, farrayPtr=maskptr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
maskptr = 0
call ESMF_GridAddCoord(make_grid, staggerloc=ESMF_STAGGERLOC_CENTER, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridAddCoord(make_grid, staggerloc=ESMF_STAGGERLOC_CORNER, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! global indexing
! atm grid is not decomposed in the y direction
!startx = lpet*atm_nx/npet*atm_dx
startx = atm_sx
starty = atm_sy
! compute coord
! X center
call ESMF_GridGetCoord(make_grid, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
coordDim=1, farrayPtr=coordX, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i = lbound(coordX,1), ubound(coordX,1)
coordX(i) = startx + atm_dx/2. + (i-1)*atm_dx
if(present(maskvalue)) then
if((coordX(i) .ge. msx .and. coordX(i) .le. mex)) then
maskptr(i,:) = maskvalue
else
maskptr(i,:) = 0
endif
endif
enddo
!print *, 'coordX: ', coordX
! X corner
call ESMF_GridGetCoord(make_grid, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
coordDim=1, farrayPtr=coordX, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i = lbound(coordX,1), ubound(coordX,1)
coordX(i) = startx + (i-1)*atm_dx
enddo
!print *, 'startx: ', startx, lbound(coordX, 1), 'coordX: ', coordX
! Y center
call ESMF_GridGetCoord(make_grid, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
coordDim=2, farrayPtr=coordY, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i = lbound(coordY,1), ubound(coordY,1)
coordY(i) = starty + atm_dy/2. + (i-1)*atm_dy
if(present(maskvalue)) then
if(.not. (coordY(i) .ge. msy .and. coordY(i) .le. mey)) then
maskptr(:,i) = 0
endif
endif
enddo
! Y corner
call ESMF_GridGetCoord(make_grid, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
coordDim=2, farrayPtr=coordY, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i = lbound(coordY,1), ubound(coordY,1)
coordY(i) = starty + (i-1)*atm_dy
enddo
if(present(field)) then
field = ESMF_FieldCreate(make_grid, typekind=ESMF_TYPEKIND_R8, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
endif
if(present(rc)) rc = ESMF_SUCCESS
end function make_grid
function make_grid_sph(atm_nx, atm_ny, atm_dx, atm_dy, atm_sx, atm_sy, &
msx, msy, mex, mey, maskvalue, tag, scheme, rc)
! return value
type(ESMF_Grid) :: make_grid_sph
! arguments
integer, intent(in) :: atm_nx, atm_ny
real(ESMF_KIND_R4), intent(in) :: atm_dx, atm_dy
real(ESMF_KIND_R4), intent(in) :: atm_sx, atm_sy
real(ESMF_KIND_R4), intent(in), optional :: msx, msy, mex, mey
integer(ESMF_KIND_I4), intent(in), optional :: maskvalue
character(len=*), intent(in), optional :: tag
integer, intent(in) , optional :: scheme
integer, intent(out), optional :: rc
! local variables
integer :: localrc, i, j
real(ESMF_KIND_R8), pointer :: coordX(:,:), coordY(:,:)
integer(ESMF_KIND_I4), pointer :: maskptr(:,:)
real(ESMF_KIND_R8) :: startx, starty
integer :: l_scheme
l_scheme = ESMF_REGRID_SCHEME_REGION3D
if(present(scheme)) l_scheme = scheme
if(l_scheme == ESMF_REGRID_SCHEME_FULL3D) then
make_grid_sph = ESMF_GridCreate1PeriDim(maxIndex=(/atm_nx, atm_ny/), &
coordSys=ESMF_COORDSYS_SPH_DEG, &
indexflag=ESMF_INDEX_GLOBAL, &
!gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/0,1/), &
!regDecomp=(/npet, 1/), &
rc=localrc)
else
make_grid_sph = ESMF_GridCreateNoPeriDim(maxIndex=(/atm_nx, atm_ny/), &
indexflag=ESMF_INDEX_GLOBAL, &
gridEdgeLWidth=(/0,0/), gridEdgeUWidth=(/1,1/), &
!regDecomp=(/npet, 1/), &
rc=localrc)
endif
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridAddItem(make_grid_sph, staggerloc=ESMF_STAGGERLOC_CENTER, &
itemflag=ESMF_GRIDITEM_MASK, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetItem(make_grid_sph, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
itemflag=ESMF_GRIDITEM_MASK, farrayPtr=maskptr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
maskptr = 0
call ESMF_GridAddCoord(make_grid_sph, staggerloc=ESMF_STAGGERLOC_CENTER, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_GridAddCoord(make_grid_sph, staggerloc=ESMF_STAGGERLOC_CORNER, &
rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! global indexing
! atm grid is not decomposed in the y direction
!startx = lpet*atm_nx/npet*atm_dx
startx = atm_sx
starty = atm_sy
! compute coord
! X center
call ESMF_GridGetCoord(make_grid_sph, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
coordDim=1, farrayPtr=coordX, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Y center
call ESMF_GridGetCoord(make_grid_sph, localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
coordDim=2, farrayPtr=coordY, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i = lbound(coordX,1), ubound(coordX,1)
do j = lbound(coordY, 2), ubound(coordY, 2)
coordX(i,j) = startx + atm_dx/2. + (i-1)*atm_dx
coordY(i,j) = starty + atm_dy/2. + (j-1)*atm_dy
enddo
enddo
if(present(maskvalue) .and. present(msx) .and. present(mex) .and. &
present(msy) .and. present(mey)) then
do i = lbound(coordX,1), ubound(coordX,1)
do j = lbound(coordY, 2), ubound(coordY, 2)
if((coordX(i,j) .ge. msx .and. coordX(i,j) .le. mex) .and. &
(coordY(i,j) .ge. msy .and. coordY(i,j) .le. mey) ) then
maskptr(i,j) = maskvalue
endif
enddo
enddo
endif
!print *, 'startx: ', startx, lbound(coordX, 1), ubound(coordX, 1), 'coordX: ', coordX(:,1)
! X corner
call ESMF_GridGetCoord(make_grid_sph, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
coordDim=1, farrayPtr=coordX, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Y corner
call ESMF_GridGetCoord(make_grid_sph, localDE=0, staggerLoc=ESMF_STAGGERLOC_CORNER, &
coordDim=2, farrayPtr=coordY, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
do i = lbound(coordX,1), ubound(coordX,1)
do j = lbound(coordX, 2), ubound(coordX, 2)
coordX(i,j) = startx + (i-1)*atm_dx
enddo
enddo
do i = lbound(coordY,1), ubound(coordY,1)
do j = lbound(coordY, 2), ubound(coordY, 2)
coordY(i,j) = starty + (j-1)*atm_dy
enddo
enddo
if(present(rc)) rc = ESMF_SUCCESS
end function make_grid_sph
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(:,:)
integer(ESMF_KIND_I4), pointer :: src_mask(:,:), dst_mask(:,:)
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*i
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
! the following no longer holds with masking, only check xgrid flux with dst field flux
!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
! mask
call ESMF_GridGetItem(dstGrid(i), localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
itemflag=ESMF_GRIDITEM_MASK, farrayPtr=dst_mask, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call compute_flux2D(vm, dst, dst_area, dst_frac, dst_frac2, allsrcsum, dst_mask, 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
!----------------------------------------------------
subroutine compute_flux1D(vm, flux_density, area, fraction, allsum, rc)
type(ESMF_VM), intent(in) :: vm
real(ESMF_KIND_R8), pointer :: flux_density(:)
real(ESMF_KIND_R8), pointer :: area(:)
real(ESMF_KIND_R8), pointer :: fraction(:)
real(ESMF_KIND_R8), intent(out) :: allsum(3)
integer, intent(out), optional :: rc
real(ESMF_KIND_R8) :: sum(3)
integer :: i,j, localrc
if(present(rc)) rc = ESMF_SUCCESS
sum = 0.
do i = lbound(flux_density, 1), ubound(flux_density, 1)
sum(1) = sum(1) + flux_density(i)*area(i)*fraction(i)
sum(2) = sum(2) + area(i)*fraction(i)
sum(3) = sum(3) + area(i)
enddo
call ESMF_VMAllReduce(vm, sum, allsum, 3, ESMF_REDUCE_SUM, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
end subroutine compute_flux1D
subroutine compute_flux2D(vm, flux_density, area, fraction, fraction2, allsum, maskptr, dstflux, rc)
type(ESMF_VM), intent(in) :: vm
real(ESMF_KIND_R8), pointer :: flux_density(:,:)
real(ESMF_KIND_R8), pointer :: area(:,:)
real(ESMF_KIND_R8), pointer :: fraction(:,:)
real(ESMF_KIND_R8), pointer :: fraction2(:,:)
real(ESMF_KIND_R8), intent(out) :: allsum(3)
integer(ESMF_KIND_I4), pointer :: maskptr(:,:)
logical, intent(in), optional :: dstflux
integer, intent(out), optional :: rc
real(ESMF_KIND_R8) :: sum(3)
integer :: i,j, localrc, npet, lpet, mask
logical :: l_dstflux
if(present(rc)) rc = ESMF_SUCCESS
l_dstflux = .false.
if(present(dstflux)) l_dstflux = dstflux
call ESMF_VMGet(vm, petCount=npet, localPet=lpet, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
!if(lpet == 0) write(*, '(A, 4I5)') 'compute_flux2D bounds: ', &
! lbound(flux_density, 1), ubound(flux_density, 1), &
! lbound(flux_density, 2), ubound(flux_density, 2)
sum = 0.
do i = lbound(flux_density, 1), ubound(flux_density, 1)
do j = lbound(flux_density, 2), ubound(flux_density, 2)
! non-masked cell has a masking value 0
if(maskptr(i,j) == 0) then
mask = 1
else
mask = 0
endif
if(l_dstflux) then
sum(1) = sum(1) + flux_density(i,j)*area(i,j)*fraction2(i,j)*mask
else
sum(1) = sum(1) + flux_density(i,j)*area(i,j)*fraction(i,j)*fraction2(i,j)*mask
endif
sum(2) = sum(2) + area(i,j)*fraction(i,j)
sum(3) = sum(3) + area(i,j)
enddo
enddo
call ESMF_VMAllReduce(vm, sum, allsum, 3, ESMF_REDUCE_SUM, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
end subroutine compute_flux2D
subroutine flux_exchange_sph(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_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(:)
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(:,:)
integer(ESMF_KIND_I4), pointer :: src_mask(:,:), dst_mask(:,:)
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
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=sideAGC
ndst=sideBGC
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, sideAGrid=srcGrid, sideBGrid=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), typekind=ESMF_TYPEKIND_R8, 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, ndst
dstField(i) = ESMF_FieldCreate(dstGrid(i), typekind=ESMF_TYPEKIND_R8, 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
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
! mask
call ESMF_GridGetItem(dstGrid(i), localDE=0, staggerLoc=ESMF_STAGGERLOC_CENTER, &
itemflag=ESMF_GRIDITEM_MASK, farrayPtr=dst_mask, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call compute_flux2D(vm, dst, dst_area, dst_frac, dst_frac2, allsrcsum, dst_mask, 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))/exf_tarea .gt. 1.e-5) .or. &
(abs(exf_tflux - allsrcsum(1))/exf_tflux .gt. 1.e-5)) 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, dst_mask, 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)/exf_tflux .gt. 1.e-5)) 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)
call ESMF_XGridDestroy(xgrid,rc=localrc)
if(present(rc)) rc = ESMF_SUCCESS
end subroutine flux_exchange_sph
end program ESMF_XGridMaskingUTest