test_cart_bilinear_xgrid Subroutine

subroutine test_cart_bilinear_xgrid(rc)

Arguments

Type IntentOptional Attributes Name
integer, intent(out) :: rc

Source Code

 subroutine test_cart_bilinear_xgrid(rc)
  integer, intent(out)  :: rc
  logical :: correct
  integer :: localrc
  type(ESMF_Grid) :: aGrid
  type(ESMF_Grid) :: bGrid
  type(ESMF_XGrid) :: dstXGrid
  type(ESMF_Mesh) :: xgridMesh
  type(ESMF_Field) :: srcField
  type(ESMF_Field) :: dstField
  type(ESMF_Field) :: xdstField
  type(ESMF_Field) :: errField
  type(ESMF_Array) :: dstArray
  type(ESMF_Array) :: errArray
  type(ESMF_Array) :: srcArray
  type(ESMF_RouteHandle) :: routeHandle
  type(ESMF_VM) :: vm
  real(ESMF_KIND_R8), pointer :: farrayPtrXC(:,:)
  real(ESMF_KIND_R8), pointer :: farrayPtrYC(:,:)
  real(ESMF_KIND_R8), pointer :: farrayPtr(:,:), farrayPtr2(:,:)
  real(ESMF_KIND_R8), pointer :: xfarrayPtr(:,:)
  real(ESMF_KIND_R8), pointer :: errfarrayPtr(:)
  real(ESMF_KIND_R8), pointer :: srcFarrayPtr(:), dstFarrayPtr(:), xdstFarrayPtr(:)
  integer :: clbnd(2),cubnd(2)
  integer :: fclbnd(2),fcubnd(2)
  integer :: i1,i2,i3, index(2)
  integer :: lDE, srclocalDECount, dstlocalDECount
  real(ESMF_KIND_R8) :: coord(2)
  character(len=ESMF_MAXSTR) :: string
  integer src_nx, src_ny, dst_nx, dst_ny
  real(ESMF_KIND_R8) :: lon, lat, theta, phi, DEG2RAD, relErr
  real(ESMF_KIND_R8) :: coords(2)
  real(ESMF_KIND_R8) :: x,y,z
  integer :: localPet, petCount, sdim
  integer :: numOwnedElems
  real(ESMF_KIND_R8), pointer :: ownedElemCoords(:)
  ! result code
  integer :: finalrc

  ! init success flag
  correct=.true.

  rc=ESMF_SUCCESS

! XMRKX !


  ! 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 = 20
  src_ny = 20

  dst_nx = 11
  dst_ny = 11

! XMRKX !
  ! degree to rad conversion
  DEG2RAD = 3.141592653589793_ESMF_KIND_R8/180.0_ESMF_KIND_R8

  ! Create A side Grid
  aGrid=ESMF_GridCreateNoPeriDimUfrm(maxIndex=(/src_nx,src_ny/), &
       minCornerCoord=(/0.0_ESMF_KIND_R8,0.0_ESMF_KIND_R8/), &
       maxCornerCoord=(/10.0_ESMF_KIND_R8,10.0_ESMF_KIND_R8/), &
       coordSys = ESMF_COORDSYS_CART, &
       staggerLocList=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER/), &
       rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif


  ! Create B side Grid
  bGrid=ESMF_GridCreateNoPeriDimUfrm(maxIndex=(/dst_nx,dst_ny/), &
       minCornerCoord=(/0.0_ESMF_KIND_R8,0.0_ESMF_KIND_R8/), &
       maxCornerCoord=(/10.0_ESMF_KIND_R8,10.0_ESMF_KIND_R8/), &
       coordSys = ESMF_COORDSYS_CART, &
       staggerLocList=(/ESMF_STAGGERLOC_CENTER, ESMF_STAGGERLOC_CORNER/), &
       rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif


  ! Create XGrid
  dstXGrid=ESMF_XGridCreate(sideAGrid=(/aGrid/),sideBGrid=(/bGrid/), &
       storeOverlay=.true., rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif


  ! Create source/destination fields
   srcField = ESMF_FieldCreate(aGrid, ESMF_TYPEKIND_R8, &
                         staggerloc=ESMF_STAGGERLOC_CENTER, name="source", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif


   dstField = ESMF_FieldCreate(dstXGrid, ESMF_TYPEKIND_R8, &
                               name="dest", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif

  xdstField = ESMF_FieldCreate(dstXGrid, ESMF_TYPEKIND_R8, &
                               name="xdest", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif

  errField = ESMF_FieldCreate(dstXGrid, ESMF_TYPEKIND_R8, &
                              name="errdest", 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(errField, array=errArray, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif


  ! srcArray
  call ESMF_FieldGet(srcField, array=srcArray, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif


  ! Get number of local DEs
  call ESMF_GridGet(aGrid, localDECount=srclocalDECount, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif


  ! Construct Src Grid
  ! (Get memory and set coords for src)
  do lDE=0,srclocalDECount-1

     ! get src pointer
     call ESMF_FieldGet(srcField, lDE, farrayPtr, &
          computationalLBound=clbnd, computationalUBound=cubnd, rc=localrc)
     if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
     endif


     !! set coords, interpolated function
     do i1=clbnd(1),cubnd(1)
     do i2=clbnd(2),cubnd(2)

        ! Get coords
        call ESMF_GridGetCoord(aGrid, staggerloc=ESMF_STAGGERLOC_CENTER, &
             localDE=lDE, index=(/i1,i2/), coord=coords, rc=localrc)
        if (localrc /=ESMF_SUCCESS) then
           rc=ESMF_FAILURE
           return
        endif

        ! set exact src data
        farrayPtr(i1,i2) = 1.0 + coords(1) + coords(2)
     enddo
     enddo

  enddo    ! lDE


  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Destination
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! Init destination field to 0.0
  ! Should only be 1 localDE
  call ESMF_FieldGet(dstField, 0, dstFarrayPtr,  rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
  endif
   
   ! Init destination field to 0.0
   dstFarrayPtr=0.0


   ! Init exact destination field
   ! Should only be 1 localDE
   call ESMF_FieldGet(xdstField, 0, xdstFarrayPtr,  rc=localrc)
   if (localrc /=ESMF_SUCCESS) then
      rc=ESMF_FAILURE
      return
   endif

  ! Get Mesh
  call ESMF_XGridGet(dstXGrid, mesh=xgridMesh, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
  endif

  ! Get number of points in destination mesh
  call ESMF_MeshGet(xgridMesh, &
       numOwnedElements=numOwnedElems, &
       spatialDim=sdim, &
       rc=localrc)
  if (ESMF_LogFoundError(localrc, &
       ESMF_ERR_PASSTHRU, &
       ESMF_CONTEXT, rcToReturn=rc)) return

   ! Allocate space for coordinates
   allocate(ownedElemCoords(sdim*numOwnedElems))


   ! Set exact destination field
   call ESMF_MeshGet(xgridMesh, ownedElemCoords=ownedElemCoords, &
       rc=localrc)
   if (ESMF_LogFoundError(localrc, &
        ESMF_ERR_PASSTHRU, &
        ESMF_CONTEXT, rcToReturn=rc)) return


   ! loop through and set xfield
   do i1=1,numOwnedElems

      ! Get coords
     coords(1)=ownedElemCoords(2*i1-1)
     coords(2)=ownedElemCoords(2*i1)

     ! Set exact dest function
     xdstFarrayPtr(i1) = 1.0 + coords(1) + coords(2)
   enddo


   ! Deallocate space for coordinates
   deallocate(ownedElemCoords)


#if 0
  call ESMF_GridWriteVTK(aGrid,staggerloc=ESMF_STAGGERLOC_CENTER, &
       filename="aGrid", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
  endif
#endif

  !!! Regrid forward from the A grid to the B grid
  ! Regrid store
  call ESMF_FieldRegridStore( &
          srcField, &
          dstField=dstField, &
          routeHandle=routeHandle, &
          regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
          rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
      rc=ESMF_FAILURE
      return
   endif

  ! 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

  ! get dst Field
  call ESMF_FieldGet(dstField, 0, dstFarrayPtr, computationalLBound=clbnd(1:1), &
                             computationalUBound=cubnd(1:1),  rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
  endif

  ! get exact destination Field
  call ESMF_FieldGet(xdstField, 0, xdstFarrayPtr,  rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
  endif

  ! Get error field
  call ESMF_FieldGet(errField, 0, errfarrayPtr,  rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
  endif


  ! destination grid
  !! check relative error
  do i1=clbnd(1),cubnd(1)

     if (xdstFarrayPtr(i1) .ne. 0.0) then
        relErr=ABS(dstFarrayPtr(i1) - xdstFarrayPtr(i1))/ABS(xdstFarrayPtr(i1))
     else
        relErr=ABS(dstFarrayPtr(i1) - xdstFarrayPtr(i1))
     endif

     ! if working everything should be close to exact answer
     if (relErr .gt. 0.005) then
        correct=.false.
        ! write(*,*) "relErr=",relErr,dstFarrayPtr(i1),xdstFarrayPtr(i1)
     endif

     ! put in error field
     errfarrayPtr(i1)=relErr

  enddo



#if 0
  call ESMF_GridWriteVTK(aGrid,staggerloc=ESMF_STAGGERLOC_CENTER, &
       filename="srcGrid", array1=srcArray, &
       rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
  endif

  call ESMF_MeshWriteVTK(xgridMesh, &
       filename="dstXGridMesh", elemarray1=dstArray, elemarray2=errArray, &
       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 xgrid
  call ESMF_XGridDestroy(dstXGrid, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
      rc=ESMF_FAILURE
      return
   endif

  ! Free the grids
  call ESMF_GridDestroy(aGrid, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
      rc=ESMF_FAILURE
      return
   endif


  call ESMF_GridDestroy(bGrid, 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_cart_bilinear_xgrid