subroutine ESMF_FieldRegridGetArea(areaField, rc)
!
! !RETURN VALUE:
!
! !ARGUMENTS:
type(ESMF_Field), intent(inout) :: areaField
integer, intent(out), optional :: rc
!
! !DESCRIPTION:
! This subroutine gets the area of the cells used for conservative interpolation for the grid object
! associated with {\tt areaField} and puts them into {\tt areaField}. If created on a 2D Grid, it must
! be built on the {\tt ESMF\_STAGGERLOC\_CENTER} stagger location.
! If created on a 3D Grid, it must be built on the {\tt ESMF\_STAGGERLOC\_CENTER\_VCENTER} stagger
! location. If created on a Mesh, it must be built on the {\tt ESMF\_MESHLOC\_ELEMENT} mesh location.
!
! If the user has set the area in the Grid, Mesh, or XGrid under {\tt areaField}, then that's the area that's
! returned in the units that the user set it in. If the user hasn't set the area, then the area is
! calculated and returned. If the Grid, Mesh, or XGrid is on the surface of a sphere, then the calculated area is in
! units of square radians. If the Grid, Mesh, or XGrid is
! Cartesian, then the calculated area is in square units of whatever unit the coordinates are in.
!
! The arguments are:
! \begin{description}
! \item [areaField]
! The Field to put the area values in.
! \item [{[rc]}]
! Return code; equals {\tt ESMF\_SUCCESS} if there are no errors.
! \end{description}
!
!EOP
integer :: localrc
type(ESMF_GeomType_Flag) :: geomtype
type(ESMF_Grid) :: Grid
type(ESMF_Array) :: Array
type(ESMF_Mesh) :: Mesh
type(ESMF_XGrid) :: xgrid
type(ESMF_StaggerLoc) :: staggerLoc, staggerLocG2M
type(ESMF_MeshLoc) :: meshloc
real(ESMF_KIND_R8), pointer :: areaFptr(:)
integer :: gridDimCount, localDECount, tileCount
type(ESMF_TypeKind_Flag) :: typekind
logical :: isLatLonDeg
integer :: isSphere
! Initialize return code; assume failure until success is certain
localrc = ESMF_SUCCESS
if (present(rc)) rc = ESMF_RC_NOT_IMPL
! Now we go through the painful process of extracting the data members
! that we need.
call ESMF_FieldGet(areaField, typekind=typekind, geomtype=geomtype, array=Array, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Check typekind
if (typekind .ne. ESMF_TYPEKIND_R8) then
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_BAD, &
msg="Area calculation is only supported for Fields of typekind=ESMF_TYPEKIND_R8", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
! If grids, then convert to a mesh to do the regridding
if (geomtype .eq. ESMF_GEOMTYPE_GRID) then
call ESMF_FieldGet(areaField, grid=Grid, &
staggerloc=staggerLoc, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Only Center stagger is supported right now until we figure out what the
! control volume for the others should be
if (staggerloc .ne. ESMF_STAGGERLOC_CENTER) then
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_BAD, &
msg="Can't currently calculate area on a stagger other then center", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
! Create the mesh from corner stagger to better represent the
! control volumes
call ESMF_GridGet(grid=Grid, tileCount=tileCount, &
dimCount=gridDimCount, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
if (gridDimCount .eq. 2) then
staggerlocG2M=ESMF_STAGGERLOC_CORNER
else if (gridDimCount .eq. 3) then
staggerlocG2M=ESMF_STAGGERLOC_CORNER_VFACE
else
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_BAD, &
msg="Can currently only do conservative regridding on 2D or 3D grids", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
! check grid
call checkGrid(Grid, staggerlocG2M, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Convert Grid to Mesh
if (tileCount .eq. 1) then
Mesh = ESMF_GridToMesh(Grid, staggerlocG2M, isSphere, isLatLonDeg, &
regridConserve=ESMF_REGRID_CONSERVE_ON, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
else
Mesh = ESMF_GridToMeshCell(Grid, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
endif
! call into the Regrid GetareaField interface
call ESMF_RegridGetArea(Grid, Mesh, Array, staggerlocG2M, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get rid of Mesh
call ESMF_MeshDestroy(Mesh, noGarbage=.true., rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
else if (geomtype .eq. ESMF_GEOMTYPE_MESH) then
call ESMF_FieldGet(areaField, mesh=Mesh, meshloc=meshloc, &
localDECount=localDECount, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
if (meshloc .ne. ESMF_MESHLOC_ELEMENT) then
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_BAD, &
msg="Can't currently calculate area on a mesh location other than elements", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
! Don't need to do anything if there are no DEs
if (localDECount < 1) then
if(present(rc)) rc = ESMF_SUCCESS
return
endif
! Get pointer to field data
! Right now Mesh will only have one DE per PET
call ESMF_FieldGet(areaField, localDE=0, farrayPtr=areaFptr, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get Area
call ESMF_MeshGetElemArea(mesh, areaList=areaFptr, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
else if (geomtype .eq. ESMF_GEOMTYPE_XGRID) then
! Get Field info
call ESMF_FieldGet(areaField, xgrid=xgrid, &
localDECount=localDECount, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Don't need to do anything if there are no DEs
if (localDECount < 1) then
if(present(rc)) rc = ESMF_SUCCESS
return
endif
! Only support 1 localDE right now
if (localDECount .ne. 1) then
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_BAD, &
msg="Getting areas for XGrid currently only supported for Fields with <= 1 local DE on each PET.", &
ESMF_CONTEXT, rcToReturn=rc)
return
endif
! Get pointer to field data
! Right now only support 1 DE per PET
call ESMF_FieldGet(areaField, localDE=0, farrayPtr=areaFptr, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Get area from XGrid
! (The ESMF_XGridGet() call checks that the size of the array matches,
! so I'm not doing it before calling in.)
call ESMF_XGridGet(xgrid, area=areaFptr, rc=localrc)
if (ESMF_LogFoundError(localrc, ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
else if (geomtype .eq. ESMF_GEOMTYPE_LOCSTREAM) then
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_BAD, &
msg="Can't get areas for a Field built on a LocStream.", &
ESMF_CONTEXT, rcToReturn=rc)
return
else
call ESMF_LogSetError(rcToCheck=ESMF_RC_ARG_BAD, &
msg="Unrecognized geometry type.", &
ESMF_CONTEXT, rcToReturn=rc)
endif
! Return success
if(present(rc)) rc = ESMF_SUCCESS
end subroutine ESMF_FieldRegridGetArea