subroutine test_regridGridToGridSph3D(rc)
integer, intent(out) :: rc
logical :: correct
integer :: localrc
type(ESMF_Grid) :: srcGrid
type(ESMF_Grid) :: dstGrid
type(ESMF_Field) :: srcField
type(ESMF_Field) :: dstField
type(ESMF_Field) :: xdstField
type(ESMF_Array) :: dstArray
type(ESMF_Array) :: xdstArray
type(ESMF_Array) :: lonArrayA
type(ESMF_Array) :: srcArrayA
type(ESMF_RouteHandle) :: routeHandle
type(ESMF_ArraySpec) :: arrayspec
type(ESMF_VM) :: vm
real(ESMF_KIND_R8), pointer :: farrayPtrXC(:,:,:), farrayPtr1D(:)
real(ESMF_KIND_R8), pointer :: farrayPtrYC(:,:,:),farrayPtrZC(:,:,:)
real(ESMF_KIND_R8), pointer :: farrayPtr(:,:,:),farrayPtr2(:,:,:)
real(ESMF_KIND_R8), pointer :: xfarrayPtr(:,:,:)
integer :: clbnd(3),cubnd(3)
integer :: fclbnd(3),fcubnd(3)
integer :: i1,i2,i3, index(3)
integer :: lDE, localDECount
real(ESMF_KIND_R8) :: coord(3)
character(len=ESMF_MAXSTR) :: string
integer dst_nx,dst_ny,dst_nz
integer src_nx,src_ny,src_nz
integer num_arrays
real(ESMF_KIND_R8) :: dst_dx,dst_dy
real(ESMF_KIND_R8) :: dst_minr, dst_maxr
real(ESMF_KIND_R8) :: src_dx,src_dy
real(ESMF_KIND_R8) :: src_minr, src_maxr
real(ESMF_KIND_R8) :: dst_minx,dst_miny,dst_minz
real(ESMF_KIND_R8) :: dst_maxx,dst_maxy,dst_maxz
real(ESMF_KIND_R8) :: x,y,z
real(ESMF_KIND_R8) :: lat_rad, lon_rad, r
real(ESMF_KIND_R8) :: err, relErr, maxRelErr
real(ESMF_KIND_R8) :: totRelErr
integer :: numRelErr
integer :: spherical_grid
integer :: localPet, petCount
real(ESMF_KIND_R8) :: beg_time, end_time
! result code
integer :: finalrc
! degree to rad conversion
real(ESMF_KIND_R8),parameter :: &
DEG2RAD = 3.141592653589793_ESMF_KIND_R8/180.0_ESMF_KIND_R8
! init success flag
correct=.true.
rc=ESMF_SUCCESS
! Init Grid resolutions
src_nx=60
src_ny=60
src_nz=14
src_minr=0.9
src_maxr=2.1
dst_nx=50
dst_ny=50
dst_nz=11
dst_minr=1.0
dst_maxr=2.0
! 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
! Setup Src Grid
srcGrid=ESMF_GridCreate1PeriDim(minIndex=(/1,1,1/),maxIndex=(/src_nx,src_ny,src_nz/), &
coordSys=ESMF_COORDSYS_SPH_DEG, regDecomp=(/petCount,1,1/), indexflag=ESMF_INDEX_GLOBAL, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Add center stagger
call ESMF_GridAddCoord(srcGrid, staggerloc=ESMF_STAGGERLOC_CENTER_VCENTER, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source field
call ESMF_ArraySpecSet(arrayspec, 3, ESMF_TYPEKIND_R8, rc=rc)
srcField = ESMF_FieldCreate(srcGrid, arrayspec, &
name="source", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Src grid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Init some stuff for convenience
src_dx = 360./src_nx
src_dy = 180./src_ny
! Get number of local DEs
call ESMF_GridGet(srcGrid, localDECount=localDECount, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get memory and set coords for dst
do lDE=0,localDECount-1
!! get coords
call ESMF_GridGetCoord(srcGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=1, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrXC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(srcGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=2, &
farrayPtr=farrayPtrYC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(srcGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=3, &
farrayPtr=farrayPtrZC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(srcField, lDE, farrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!! set coords
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
do i3=clbnd(3),cubnd(3)
! Set coordinates
farrayPtrXC(i1,i2,i3) = REAL(i1-1)*src_dx
farrayPtrYC(i1,i2,i3) = -90. + (REAL(i2-1)*src_dy + 0.5*src_dy)
farrayPtrZC(i1,i2,i3) = ((src_maxr-src_minr)*REAL(i3-1)/REAL(src_nz-1))+src_minr
! Compute Lat/Lon in rad
lon_rad = DEG2RAD*farrayPtrXC(i1,i2,i3)
lat_rad = DEG2RAD*(90.-farrayPtrYC(i1,i2,i3))
r=farrayPtrZC(i1,i2,i3)
! initialize source field
farrayPtr(i1,i2,i3)= r*5.0+2.0+ cos(lon_rad)**2.*cos(2.*lat_rad)
!farrayPtr(i1,i2,i3)= 2.0+r*5.0+ cos(lat_rad)**2.*cos(2.*lon_rad)
!farrayPtr(i1,i2,i3)= 1.0
enddo
enddo
enddo
enddo ! lDE
! Setup Dst Grid
dstGrid=ESMF_GridCreate1PeriDim(minIndex=(/1,1,1/),maxIndex=(/dst_nx,dst_ny,dst_nz/), &
coordSys=ESMF_COORDSYS_SPH_DEG, regDecomp=(/1,petCount,1/), indexflag=ESMF_INDEX_GLOBAL, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Create source/destination fields
call ESMF_ArraySpecSet(arrayspec, 3, ESMF_TYPEKIND_R8, rc=rc)
dstField = ESMF_FieldCreate(dstGrid, arrayspec, &
staggerloc=ESMF_STAGGERLOC_CENTER_VCENTER, name="dest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Exact destination field
xdstField = ESMF_FieldCreate(dstGrid, arrayspec, &
staggerloc=ESMF_STAGGERLOC_CENTER_VCENTER, name="dest", rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridAddCoord(dstGrid, staggerloc=ESMF_STAGGERLOC_CENTER_VCENTER, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get arrays
! dstArray
call ESMF_FieldGet(dstField, array=dstArray, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(xdstField, array=xdstArray, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! srcArrayA
call ESMF_FieldGet(srcField, array=srcArrayA, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Destination grid
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Init some stuff for convenience
dst_dx = 360./dst_nx
dst_dy = 180./dst_ny
! Get number of local DEs
call ESMF_GridGet(dstGrid, localDECount=localDECount, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get memory and set coords for dst
do lDE=0,localDECount-1
!! get coords
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=1, &
computationalLBound=clbnd, computationalUBound=cubnd, farrayPtr=farrayPtrXC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=2, &
farrayPtr=farrayPtrYC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER_VCENTER, coordDim=3, &
farrayPtr=farrayPtrZC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(dstField, lDE, farrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(xdstField, lDE, xfarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!! set coords
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
do i3=clbnd(3),cubnd(3)
! Set coordinates
farrayPtrXC(i1,i2,i3) = REAL(i1-1)*dst_dx
farrayPtrYC(i1,i2,i3) = -90. + (REAL(i2-1)*dst_dy + 0.5*dst_dy)
farrayPtrZC(i1,i2,i3) = ((dst_maxr-dst_minr)*REAL(i3-1)/REAL(dst_nz-1))+dst_minr
! Compute Lat/Lon in rad
lon_rad = DEG2RAD*farrayPtrXC(i1,i2,i3)
lat_rad = DEG2RAD*(90.-farrayPtrYC(i1,i2,i3))
r=farrayPtrZC(i1,i2,i3)
! initialize exact destination field
!xfarrayPtr(i1,i2,i3)=r*(2. + cos(lon_rad)**2.*cos(2.*lat_rad))
!xfarrayPtr(i1,i2,i3)=r*5.0+2. + cos(lon_rad)**2.*cos(2.*lat_rad)
!xfarrayPtr(i1,i2,i3)= 1.0+r*0.01+ cos(lat_rad)**2.*cos(2.*lon_rad)
xfarrayPtr(i1,i2,i3)= r*5.0+2.0+ cos(lon_rad)**2.*cos(2.*lat_rad)
!xfarrayPtr(i1,i2,i3)= 1.0
!xfarrayPtr(i1,i2,i3)=r
! Initialize destination field
farrayPtr(i1,i2,i3)=0.0
enddo
enddo
enddo
enddo ! lDE
! Get start time
! call ESMF_VMWtime(beg_time)
!!! Regrid forward from the A grid to the B grid
! Regrid store
call ESMF_FieldRegridStore( &
srcField, &
dstField=dstField, &
routeHandle=routeHandle, &
regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
lineType=ESMF_LINETYPE_GREAT_CIRCLE, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Get end time
!call ESMF_VMWtime(end_time)
! output info
! write(*,*) "Src dims=",src_nx,src_ny
! write(*,*) "Dst dims=",dst_nx,dst_ny
! write(*,*) "Store time = ",end_time-beg_time
! Do regrid
call ESMF_FieldRegrid(srcField, dstField, 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
! Init max
maxRelErr=0.0
totRelErr=0.0
numRelErr=0
! Check error
do lDE=0,localDECount-1
!! get coords
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=1, &
farrayPtr=farrayPtrXC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=2, &
farrayPtr=farrayPtrYC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridGetCoord(dstGrid, localDE=lDE, staggerLoc=ESMF_STAGGERLOC_CENTER, coordDim=3, &
farrayPtr=farrayPtrZC, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(dstField, lDE, farrayPtr, computationalLBound=clbnd, &
computationalUBound=cubnd, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldGet(xdstField, lDE, xfarrayPtr, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
!! check error
do i1=clbnd(1),cubnd(1)
do i2=clbnd(2),cubnd(2)
do i3=clbnd(3),cubnd(3)
! Calculate error
err=abs(farrayPtr(i1,i2,i3)-xfarrayPtr(i1,i2,i3))
! Calculate relative error
if (xfarrayPtr(i1,i2,i3) .ne. 0.0) then
relErr=err/xfarrayPtr(i1,i2,i3)
else
relErr=err
endif
!! Calculate max
if (relErr > maxRelErr) maxRelErr=relErr
!! Accumulate avg. info
totRelErr=totRelErr+maxRelErr
numRelErr = numRelErr+1
!! if error is too big report an error
if (relErr > 0.002) then
correct=.false.
endif
enddo
enddo
enddo
enddo ! lDE
if (.not. correct) then
write(*,*) "Test not correct. Max Rel. Error= ",maxRelErr
endif
!write(*,*) "Max Rel. Error= ",maxRelErr
!write(*,*) "Avg Rel. Error= ",totRelErr/REAL(numRelErr)
#if 0
call ESMF_GridWriteVTK(srcGrid,staggerloc=ESMF_STAGGERLOC_CENTER, &
filename="srcGrid", &
rc=localrc)
call ESMF_GridWriteVTK(dstGrid,staggerloc=ESMF_STAGGERLOC_CENTER, &
filename="dstGrid", array1=dstArray, array2=xdstArray, &
rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
#endif
! Destroy the Fields
call ESMF_FieldDestroy(srcField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_FieldDestroy(dstField, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
! Free the grids
call ESMF_GridDestroy(srcGrid, rc=localrc)
if (localrc /=ESMF_SUCCESS) then
rc=ESMF_FAILURE
return
endif
call ESMF_GridDestroy(dstGrid, 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_regridGridToGridSph3D