subroutine test_regrid180vs360(rc)
integer, intent(out) :: rc
logical :: correct
integer :: localrc
type(ESMF_Grid) :: grid360
type(ESMF_Grid) :: grid180
type(ESMF_Field) :: srcField360
type(ESMF_Field) :: dstField360
type(ESMF_Field) :: field180
type(ESMF_Field) :: errorField
type(ESMF_Array) :: array180
type(ESMF_Array) :: errorArray
type(ESMF_Array) :: lonArray360
type(ESMF_Array) :: srcArray360, dstArray360
type(ESMF_RouteHandle) :: routeHandle
type(ESMF_RouteHandle) :: routeHandle1
type(ESMF_ArraySpec) :: arrayspec
type(ESMF_VM) :: vm
real(ESMF_KIND_R8), pointer :: farrayPtrXC(:,:)
real(ESMF_KIND_R8), pointer :: farrayPtrYC(:,:)
real(ESMF_KIND_R8), pointer :: farrayPtr(:,:),farrayPtr2(:,:),errorfarrayPtr(:,:)
integer :: petMap2D(2,2,1)
integer :: clbnd(2),cubnd(2)
integer :: fclbnd(2),fcubnd(2)
integer :: i1,i2, index(2)
integer :: lDE, localDECount
real(ESMF_KIND_R8) :: coord(2)
character(len=ESMF_MAXSTR) :: string
integer src_nx, src_ny, dst_nx, dst_ny
integer num_arrays
real(ESMF_KIND_R8) :: src_dx, src_dy
real(ESMF_KIND_R8) :: dst_dx, dst_dy
real(ESMF_KIND_R8) :: ctheta, stheta
real(ESMF_KIND_R8) :: theta, d2rad, x, y, z
real(ESMF_KIND_R8) :: DEG2RAD, a, lat, lon, phi
real(ESMF_KIND_R8) :: rangle, xtmp, ytmp, ztmp
real(ESMF_KIND_R8) :: RAD2DEG
integer :: spherical_grid
integer, pointer :: larrayList(:)
integer :: localPet, petCount
integer :: srcTermProcessing, pipeLineDepth
! result code
integer :: finalrc
! init success flag
correct=.true.
rc=ESMF_SUCCESS
! get pet info
call ESMF_VMGetGlobal(vm, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
call ESMF_VMGet(vm, petCount=petCount, localPet=localpet, rc=localrc)
if (ESMF_LogFoundError(localrc, &
ESMF_ERR_PASSTHRU, &
ESMF_CONTEXT, rcToReturn=rc)) return
! Establish the resolution of the grids
src_nx = 100
src_ny = 50
dst_nx = 90
dst_ny = 40
! setup source grid
grid360=ESMF_GridCreate1PeriDim(minIndex=(/1,1/),maxIndex=(/src_nx,src_ny/),regDecomp=(/petCount,1/), &
coordSys=ESMF_COORDSYS_SPH_DEG, indexflag=ESMF_INDEX_GLOBAL, &
! poleType=(/ESMF_POLETYPE_MONOPOLE,ESMF_POLETYPE_BIPOLE/), &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! setup dest. grid
grid180=ESMF_GridCreate1PeriDim(minIndex=(/1,1/),maxIndex=(/dst_nx,dst_ny/),regDecomp=(/1,petCount/), &
coordSys=ESMF_COORDSYS_SPH_DEG, indexflag=ESMF_INDEX_GLOBAL, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source/destination fields
call ESMF_ArraySpecSet(arrayspec, 2, ESMF_TYPEKIND_R8, rc=rc)
srcField360 = ESMF_FieldCreate(grid360, arrayspec, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="source", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
dstField360 = ESMF_FieldCreate(grid360, arrayspec, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="source", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
errorField = ESMF_FieldCreate(grid360, arrayspec, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="source", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
field180 = ESMF_FieldCreate(grid180, arrayspec, &
staggerloc=ESMF_STAGGERLOC_CENTER, name="dest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Allocate coordinates
call ESMF_GridAddCoord(grid360, staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridAddCoord(grid180, staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get number of local DEs
call ESMF_GridGet(grid360, localDECount=localDECount, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get arrays
! array180
call ESMF_FieldGet(field180, array=array180, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! srcArray360
call ESMF_FieldGet(srcField360, array=srcArray360, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! dstArray360
call ESMF_FieldGet(dstField360, array=dstArray360, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! errorArray
call ESMF_FieldGet(errorField, array=errorArray, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!! get longitude array
call ESMF_GridGetCoord(grid360, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, &
array=lonArray360, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Write results to a mesh
num_arrays = 1
! Test interpolation on the sphere
! Set the source grid coordinates to be a 0 to 360 grid
src_dx = 360./src_nx
src_dy = 180./src_ny
DEG2RAD = 3.14159265/180.0
RAD2DEG = 1./DEG2RAD
! Get memory and set coords for src
do lDE=0,localDECount-1
!! get coord 1
call ESMF_GridGetCoord(grid360, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrXC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(grid360, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrYC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get src pointer
call ESMF_FieldGet(srcField360, lDE, farrayPtr, computationalLBound=fclbnd, &
computationalUBound=fcubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get src destination Field pointer
call ESMF_FieldGet(dstField360, lDE, farrayPtr2, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
if (clbnd(1) .ne. fclbnd(1)) print *, 'Error clbnd != fclbnd'
if (clbnd(2) .ne. fclbnd(2)) print *, 'Error clbnd != fclbnd'
if (cubnd(1) .ne. fcubnd(1)) print *, 'Error cubnd != fcubnd'
if (cubnd(2) .ne. fcubnd(2)) print *, 'Error cubnd != fcubnd'
!! set coords, interpolated function
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
! Set source coordinates as 0 to 360
farrayPtrXC(i1,i2) = REAL(i1-1)*src_dx
farrayPtrYC(i1,i2) = -90. + (REAL(i2-1)*src_dy + 0.5*src_dy)
lon = farrayPtrXC(i1,i2)
lat = farrayPtrYC(i1,i2)
! Set the source to be a function of the x,y,z coordinate
theta = DEG2RAD*(lon)
phi = DEG2RAD*(90.-lat)
x = cos(theta)*sin(phi)
y = sin(theta)*sin(phi)
z = cos(phi)
! set src data
! (something relatively smooth, that varies everywhere)
farrayPtr(i1,i2) = x+y+z+15.0
! initialize src destination field
farrayPtr2(i1,i2)=0.0
enddo
enddo
enddo ! lDE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Destination grid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
dst_dx = 360./dst_nx
dst_dy = 180./dst_ny
rangle = DEG2RAD*20.
! Get memory and set coords for dst
do lDE=0,localDECount-1
!! get coord 1
call ESMF_GridGetCoord(grid180, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrXC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(grid180, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrYC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(field180, lDE, farrayPtr, computationalLBound=fclbnd, &
computationalUBound=fcubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
if (clbnd(1) .ne. fclbnd(1)) print *, 'Error dst clbnd != fclbnd'
if (clbnd(2) .ne. fclbnd(2)) print *, 'Error dst clbnd != fclbnd'
if (cubnd(1) .ne. fcubnd(1)) print *, 'Error dst cubnd != fcubnd'
if (cubnd(2) .ne. fcubnd(2)) print *, 'Error dst cubnd != fcubnd'
!! set coords, interpolated function
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
! Set destination coordinates as -180 to 180
farrayPtrXC(i1,i2) = -180. + (REAL(i1-1)*dst_dx)
farrayPtrYC(i1,i2) = -90. + (REAL(i2-1)*dst_dy + 0.5*dst_dy)
! init destination mesh to 0
farrayPtr(i1,i2) = 0.
enddo
enddo
enddo ! lDE
!!! Regrid forward from the 0 to 360 grid to the -180 to 180 grid
! Regrid store
! make sure these work
srcTermProcessing=0
pipeLineDepth=localPet
call ESMF_FieldRegridStore(srcField360, dstField=field180, &
routeHandle=routeHandle, &
regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
srcTermProcessing=srcTermProcessing, &
pipeLineDepth=pipeLineDepth, &
checkFlag=.true., & ! Just add this to make sure it doesn't cause problems, even when no error
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Do regrid
call ESMF_FieldRegrid(srcField360, field180, routeHandle, &
termorderflag=ESMF_TERMORDER_FREE ,rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldRegridRelease(routeHandle, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!!!!!!!! Regrid back from the -180 to 180 grid to the 0 to 360 grid
! Regrid store
call ESMF_FieldRegridStore(field180, dstField=dstField360, &
routeHandle=routeHandle, &
regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Do regrid
call ESMF_FieldRegrid(field180, dstField360, routeHandle, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldRegridRelease(routeHandle, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Check if the values are close
do lDE=0,localDECount-1
! get src Field
call ESMF_FieldGet(srcField360, lDE, farrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get src destination Field
call ESMF_FieldGet(dstField360, lDE, farrayPtr2, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! get src Field
call ESMF_FieldGet(errorField, lDE, errorfarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!! check relative error
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
if (farrayPtr(i1,i2) .ne. 0.0) then
errorfarrayPtr(i1,i2)=ABS((farrayPtr(i1,i2) - farrayPtr2(i1,i2))/farrayPtr(i1,i2))
else
errorfarrayPtr(i1,i2)=(farrayPtr(i1,i2) - farrayPtr2(i1,i2))
endif
if (ABS(errorfarrayPtr(i1,i2)) .gt. 0.01) then
correct=.false.
endif
enddo
enddo
enddo ! lDE
! Uncomment these calls to see some actual regrid results
#if 0
spherical_grid = 1
call ESMF_MeshIO(vm, grid360, ESMF_STAGGERLOC_CENTER, &
"srcmesh", srcArray360, dstArray360, errorArray, lonArray360, rc=localrc, &
spherical=spherical_grid)
write(*,*) "LOCALRC=",localrc
call ESMF_MeshIO(vm, grid180, ESMF_STAGGERLOC_CENTER, &
"dstmesh", array180, rc=localrc, &
spherical=spherical_grid)
write(*,*) "LOCALRC=",localrc
#endif
! Destroy the Fields
call ESMF_FieldDestroy(srcField360, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(dstField360, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(errorField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(field180, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Free the grids
call ESMF_GridDestroy(grid360, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridDestroy(grid180, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! return answer based on correct flag
if (correct) then
rc=ESMF_SUCCESS
else
rc=ESMF_FAILURE
endif
end subroutine test_regrid180vs360