subroutine ESMF_RegridWeightGenCheck(weightFile, keywordEnforcer, &
checkMethod, rc)
! !ARGUMENTS:
character(len=*), intent(in) :: weightFile
type(ESMF_KeywordEnforcer), optional:: keywordEnforcer ! must use keywords below
type(ESMF_RWGCheckMethod_Flag), intent(in), optional :: checkMethod
integer, intent(out), optional :: rc
! !DESCRIPTION:
!
! The arguments are:
! \begin{description}
! \item [weightFile]
! The interpolation weight file name.
! \item [{[checkMethod]}]
! Method to use when checking sparse matrix multiplication. There are two
! options. {\tt ESMF\_CHECKMETHOD\_FIELD} (the default) uses field sparse
! matrix multiplication. {\tt ESMF\_CHECKMETHOD\_ARRAY} uses array sparse
! matrix multiplication. The underlying sparse matrix multiplication code
! does not change between field and array. The difference is in how the
! higher level objects call into the sparse matrix multiplication.
! \item [{[rc]}]
! Return code; equals {\tt ESMF\_SUCCESS} if there are no errors.
! \end{description}
!EOPI
!--------------------------------------------------------------------------
! DECLARATIONS
!--------------------------------------------------------------------------
integer :: localPet, nPet
integer :: status
type(ESMF_VM) :: vm
character(ESMF_MAXPATHLEN) :: title
real(ESMF_KIND_R8), allocatable :: factorList(:)
integer, allocatable :: factorIndexList(:,:)
real(ESMF_KIND_R8), pointer :: src_lat(:), src_lon(:), &
dst_lat(:), dst_lon(:), &
src_area(:), dst_area(:), &
src_mask(:), dst_mask(:), &
src_frac(:), dst_frac(:)
integer :: src_dim, dst_dim
integer :: i, src, dst
real(ESMF_KIND_R8), parameter :: two = 2.0
real(ESMF_KIND_R8), parameter :: d2r = 3.141592653589793238/180
real(ESMF_KIND_R8), parameter :: UNINITVAL = 422397696
real(ESMF_KIND_R8), allocatable :: FsrcArray(:)
real(ESMF_KIND_R8), allocatable :: FdstArray(:), FdstArrayX(:)
real(ESMF_KIND_R8), pointer :: FdstArrayPtr(:), farrayPtr(:)
type(ESMF_DistGrid) :: src_distgrid, dst_distgrid
type(ESMF_ArraySpec):: src_arrayspec, dst_arrayspec
type(ESMF_Array) :: srcArray, dstArray
type(ESMF_RouteHandle) :: routehandle
type(ESMF_Grid) :: srcGrid, dstGrid
type(ESMF_Field) :: srcField, dstField
real(ESMF_KIND_R8) :: maxRelError, meanRelError, lsRelError
real(ESMF_KIND_R8) :: totErrDif, twoErrDif, twoErrX
real(ESMF_KIND_R8) :: err, relError, maxneg, maxpos
integer:: numRelError
real(ESMF_KIND_R8) :: grid1min, grid1max, grid2min, grid2max
real(ESMF_KIND_R8), pointer :: grid1area(:), grid2area(:)
real(ESMF_KIND_R8), pointer :: grid1areaXX(:), grid2areaXX(:)
type(ESMF_NormType_Flag) :: normType
type(ESMF_RegridMethod_Flag) :: regridmethod
type(ESMF_RWGCheckMethod_Flag) :: checkMethodLocal
!--------------------------------------------------------------------------
! EXECUTION
!--------------------------------------------------------------------------
! Set the default for the SMM check method.
if (.not. present(checkMethod)) then
checkMethodLocal = ESMF_RWGCHECKMETHOD_FIELD
else
checkMethodLocal = checkMethod
endif
#ifdef ESMF_NETCDF
! set log to flush after every message
call ESMF_LogSet(flush=.true., rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! get all vm information
call ESMF_VMGetGlobal(vm, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! set up local pet info
call ESMF_VMGet(vm, localPet=localPet, petCount=nPet, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
!Set finalrc to success
rc = ESMF_SUCCESS
! read in the grid dimensions
call NCFileInquire(weightFile, title, normType, src_dim, dst_dim, &
regridmethod, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! only read the data on PET 0 until we get ArrayRead going...
if (localPet == 0) then
allocate(src_lat(src_dim))
allocate(src_lon(src_dim))
allocate(src_area(src_dim))
allocate(src_mask(src_dim))
allocate(src_frac(src_dim))
allocate(dst_lat(dst_dim))
allocate(dst_lon(dst_dim))
allocate(dst_area(dst_dim))
allocate(dst_mask(dst_dim))
allocate(dst_frac(dst_dim))
call GridReadCoords(weightFile, src_lat, src_lon, src_area, src_mask, src_frac, &
dst_lat, dst_lon, dst_area, dst_mask, dst_frac, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! create Fortran arrays
allocate(FsrcArray(src_dim))
allocate(FdstArray(dst_dim))
allocate(FdstArrayX(dst_dim))
! Initialize data
FsrcArray = two + cos(src_lat)**2*cos(two*src_lon)
FdstArrayX = two + cos(dst_lat)**2*cos(two*dst_lon)
FdstArray = UNINITVAL
! deallocate arrays
deallocate(src_lat)
deallocate(src_lon)
deallocate(dst_lat)
deallocate(dst_lon)
endif
! create DistGrids for the ESMF Arrays
src_distgrid = ESMF_DistGridCreate(minIndex=(/1/), &
maxIndex=(/src_dim/), rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
dst_distgrid = ESMF_DistGridCreate(minIndex=(/1/), &
maxIndex=(/dst_dim/), rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! create dummy grids for fields
srcGrid = ESMF_GridCreate(src_distgrid, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
dstGrid = ESMF_GridCreate(dst_distgrid, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! create ArraySpecs for the ESMF Arrays
call ESMF_ArraySpecSet(src_arrayspec, typekind=ESMF_TYPEKIND_R8, rank=1, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
call ESMF_ArraySpecSet(dst_arrayspec, typekind=ESMF_TYPEKIND_R8, rank=1, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! create the ESMF Arrays
srcArray = ESMF_ArrayCreate(arrayspec=src_arrayspec, distgrid=src_distgrid, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
dstArray = ESMF_ArrayCreate(arrayspec=dst_arrayspec, distgrid=dst_distgrid, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! initialize the array..
call ESMF_ArrayGet(dstArray, farrayPtr=farrayPtr, rc=rc)
farrayPtr = UNINITVAL
! Scatter the ESMF Arrays
call ESMF_ArrayScatter(srcArray, farray=FsrcArray, rootPet=0, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! create fields on the empty grid and arrays
srcField = ESMF_FieldCreate(srcGrid, srcArray, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! create fields on the empty grid and arrays
dstField = ESMF_FieldCreate(dstGrid, dstArray, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
#if 0
!-------------------- diagnostics -----------------------------------------
print *, "size of factorList = (", size(factorList,1),")"
print *, "size of factorIndexList = (", size(factorIndexList,1), ",", size( factorIndexList,2),")"
do i=1,size(factorList)
src = factorIndexList(1,i)
dst = factorIndexList(2,i)
FdstArray(dst) = FdstArray(dst) + factorList(i)*FsrcArray(src)
!print *, FdstArray(dst), " ", FsrcArray(factorIndexList(1,i)), " ", factorList(i)
!print *, factorList(i), FsrcArray(factorIndexList(1,i)), " ", FdstArrayX( factorIndexList(2,i))
enddo
#endif
if (checkMethodLocal .eq. ESMF_RWGCHECKMETHOD_FIELD) then
! Field and Grid way of doing things
call ESMF_FieldSMMStore(srcField=srcField, dstField=dstField, &
filename=weightFile, routehandle=routehandle, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! compute a Regrid from srcField to dstField
call ESMF_FieldRegrid(srcField, dstField, routehandle, &
zeroregion=ESMF_REGION_SELECT, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
call ESMF_FieldRegridRelease(routehandle, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
else if (checkMethodLocal .eq. ESMF_RWGCHECKMETHOD_ARRAY) then
! Array way of doing things
call ESMF_ArraySMMStore(srcArray=srcArray, dstArray=dstArray, &
filename=weightFile, routehandle=routehandle, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! *************************************************************
! compute a SMM from srcArray to dstArray
call ESMF_ArraySMM(srcArray, dstArray, routehandle, &
zeroregion=ESMF_REGION_SELECT, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
call ESMF_ArraySMMRelease(routehandle, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
else
if (ESMF_LogFoundError(ESMF_RC_ARG_BAD, ESMF_ERR_PASSTHRU, ESMF_CONTEXT, &
rcToReturn=rc)) return
endif
! ArrayGather the dst array
call ESMF_ArrayGather(dstArray, farray=FdstArray, rootPet=0, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! -----------------------------------------------------------------------
! ERROR ANALYSIS - serial
! -----------------------------------------------------------------------
if (localPet == 0) then
totErrDif = 0
twoErrDif = 0
twoErrX = 0
maxRelError=0.0
meanRelError=0.0
lsRelError=0.0
numRelError=0
! source error
grid1min = UNINITVAL
grid1max = 0
do i=1,src_dim
if (src_mask(i) /=0) then
if(FsrcArray(i) < grid1min) grid1min = FsrcArray(i)
if(FsrcArray(i) > grid1max) grid1max = FsrcArray(i)
endif
enddo
! destination error
grid2min = UNINITVAL
grid2max = 0
do i=1,dst_dim
! don't look in masked cells
! if frac is below .999, then a significant portion of this cell is
! missing from the weight calculation and error is misleading here
! also don't look in unitialized cells, for the regional to global cases
if (dst_mask(i) /= 0 .and. dst_frac(i) > .999 &
.and. FdstArray(i) /= UNINITVAL) then
! compute the raw pointwise error
err = abs(FdstArray(i) - FdstArrayX(i))
! Compute relative error
if (FdstArrayX(i) .ne. 0.0) then
relError=err/FdstArrayX(i)
else
relError=err
endif
! Compute the max relative error
if (relError > maxRelError) then
maxRelError = relError
endif
! relative pointwise error summation
totErrDif = totErrDif + relError
twoErrDif = twoErrDif + relError**2
twoErrX = twoErrX + FdstArrayX(i)**2
! raw pointwise error summation
!totErrDif = totErrDif + err
!twoErrDif = twoErrDif + err**2
!twoErrX = twoErrX + FdstArrayX(i)**2
! Number of points with frac > .999
numRelError=numRelError+1
! masking will screw this one up
if (FdstArray(i) < grid2min) grid2min = FdstArray(i)
if (FdstArray(i) > grid2max) grid2max = FdstArray(i)
endif
enddo
! Read factors from weights file. Factors are allocated in the
! subroutine. Only the factorList is needed for error assessment.
call ESMF_FactorRead(weightFile, factorList, factorIndexList, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
! maximum negative weight
maxneg = 0
maxneg = minval(factorList)
if (maxneg > 0) maxneg = 0
! maximum positive weight
maxpos = 0
maxpos = maxval(factorList)
deallocate(factorList, factorIndexList)
! error measures
meanRelError = totErrDif/REAL(numRelError,ESMF_KIND_R8)
lsRelError = sqrt(twoErrDif)/sqrt(twoErrX)
! area calculations for conservative regridding only
! Calculate dst area depending on norm option
if (regridmethod == ESMF_REGRIDMETHOD_CONSERVE) then
! use one of src_ or dst_frac, but NOT both!
allocate(grid1area(src_dim))
allocate(grid2area(dst_dim))
allocate(grid1areaXX(src_dim))
allocate(grid2areaXX(dst_dim))
grid1area = FsrcArray*src_area*src_frac
grid2area=0.0
if (normType == ESMF_NORMTYPE_DSTAREA) then
do i=1,dst_dim
! Only calculate dst area over region that is unmasked and initialized
if ((dst_mask(i) /= 0) .and. (FdstArray(i) /=UNINITVAL)) then
grid2area(i) = FdstArray(i)*dst_area(i)
endif
enddo
else if (normType == ESMF_NORMTYPE_FRACAREA) then
do i=1,dst_dim
! Only calculate dst area over region that is unmasked and initialized
if ((dst_mask(i) /= 0) .and. (FdstArray(i) /=UNINITVAL)) then
grid2area(i) = FdstArray(i)*dst_area(i)*dst_frac(i)
endif
enddo
endif
grid1areaXX = FsrcArray*src_area
grid2areaXX = FdstArray*dst_area*dst_frac
endif
print *, trim(weightFile), " - ", trim(title)
print *, " "
print *, "Grid 1 min: ", grid1min, " Grid 1 max: ", grid1max
print *, "Grid 2 min: ", grid2min, " Grid 2 max: ", grid2max
print *, " "
print *, "Maximum negative weight = ", maxneg
print *, "Maximum positive weight = ", maxpos
print *, " "
print *, "Mean relative error = ", meanRelError
print *, "Maximum relative error = ", maxRelError
print *, "Least squares error = ", lsRelError
print *, " "
if (regridmethod == ESMF_REGRIDMETHOD_CONSERVE) then
print *, "Grid 1 area = ", sum(grid1area)
print *, "Grid 2 area = ", sum(grid2area)
print *, "Conservation error = ", abs(sum(grid2area)-sum(grid1area))
deallocate(grid1area)
deallocate(grid2area)
deallocate(grid1areaXX)
deallocate(grid2areaXX)
else
print *, "Grid 1 area = 0.000000000000000"
print *, "Grid 2 area = 0.000000000000000"
print *, "Conservation error = 0.000000000000000"
endif
! print *, " "
! print *, "reverse fracs - Grid 1 area = ", sum(grid1areaXX)
! print *, "reverse fracs - Grid 2 area = ", sum(grid2areaXX)
! print *, "reverse - Conservation error = ", abs(sum(grid2areaXX)-sum(grid1areaXX))
deallocate(src_area)
deallocate(src_frac)
deallocate(dst_area)
deallocate(dst_frac)
deallocate(src_mask)
deallocate(dst_mask)
deallocate(FsrcArray)
deallocate(FdstArray)
deallocate(FdstArrayX)
endif
! destroy and deallocate
call ESMF_ArrayDestroy(srcArray, rc=status)
call ESMF_ArrayDestroy(dstArray, rc=status)
call ESMF_GridDestroy(srcGrid, rc=status)
call ESMF_GridDestroy(dstGrid, rc=status)
call ESMF_FieldDestroy(srcField, rc=status)
call ESMF_FieldDestroy(dstField, rc=status)
if (ESMF_LogFoundError(status, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, &
rcToReturn=rc)) return
#else
call ESMF_LogSetError(rcToCheck=ESMF_RC_LIB_NOT_PRESENT, &
msg="- ESMF_NETCDF not defined when lib was compiled", &
ESMF_CONTEXT, rcToReturn=rc)
return
#endif
end subroutine ESMF_RegridWeightGenCheck