test_regridNearestLocStreamToGrid Subroutine

subroutine test_regridNearestLocStreamToGrid(rc)

Arguments

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

Source Code

 subroutine test_regridNearestLocStreamToGrid(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) :: arrayB
  type(ESMF_Array) :: srcArrayA
  type(ESMF_RouteHandle) :: routeHandle
  type(ESMF_ArraySpec) :: arrayspec1,arrayspec2
  type(ESMF_VM) :: vm
  real(ESMF_KIND_R8), pointer :: farrayPtrXC(:,:)
  real(ESMF_KIND_R8), pointer :: farrayPtrYC(:,:)
  real(ESMF_KIND_R8), pointer :: srcPtr(:),dstPtr(:,:)
  real(ESMF_KIND_R8), pointer :: xdstPtr(:,:)
  integer :: clbnd(2),cubnd(2)
  integer :: fclbnd(2),fcubnd(2)
  integer :: i1,i2,i3, index(2)
  integer :: lDE, localDECount
  integer :: cl,cu,cc
  integer :: clb,cub,mcc,elb,eub,ec
  real(ESMF_KIND_R8) :: coord(2)
  character(len=ESMF_MAXSTR) :: string
  integer src_nx, src_ny, dst_nx, dst_ny
  integer :: numLocationsOnThisPet,idx
  type(ESMF_LocStream) :: srcLocStream
  real(ESMF_KIND_R8), pointer :: Xarray(:),Yarray(:)

  real(ESMF_KIND_R8) :: theta, phi
  real(ESMF_KIND_R8) :: src_dx, src_dy
  real(ESMF_KIND_R8) :: dst_dx, dst_dy
    ! degree to rad conversion
  real(ESMF_KIND_R8),parameter :: DEG2RAD = 3.141592653589793_ESMF_KIND_R8/180.0_ESMF_KIND_R8
  integer :: localPet, petCount

  ! result code
  integer :: finalrc

  ! init flags
  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

  src_nx = 20
  src_ny = 20
  dst_nx = 20
  dst_ny = 20

  src_dx=360.0/src_nx
  src_dy=180.0/src_ny
  dst_dx=360.0/dst_nx
  dst_dy=180.0/dst_ny


  srcLocStream=ESMF_LocStreamCreate(minIndex=1, maxIndex=src_nx*src_ny, regDecomp=petCount, &
                                    indexflag=ESMF_INDEX_GLOBAL, &
                                    coordSys=ESMF_COORDSYS_SPH_DEG, &
                                    rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    print*,'ERROR:  trouble creating locStream'
    rc=ESMF_FAILURE
    return
  endif

  call  ESMF_LocStreamGet(srcLocStream, localDECount=localDECount, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    print*,'ERROR:  trouble with locStreamGet'
    rc=ESMF_FAILURE
    return
  endif

  !-------------------------------------------------------------------
  ! Add key data (internally allocating memory).
  !-------------------------------------------------------------------
  call ESMF_LocStreamAddKey(srcLocStream,                 &
                            keyName="ESMF:Lat",           &
                            KeyTypeKind=ESMF_TYPEKIND_R8, &
                            keyUnits="degrees",           &
                            keyLongName="Latitude", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    print*,'ERROR:  trouble adding LocStream key for latitude'
    rc=ESMF_FAILURE
    return
  endif
  call ESMF_LocStreamAddKey(srcLocStream,                 &
                            keyName="ESMF:Lon",           &
                            KeyTypeKind=ESMF_TYPEKIND_R8, &
                            keyUnits="degrees",           &
                            keyLongName="Longitude", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    print*,'ERROR:  trouble adding LocStream key for longitude'
    rc=ESMF_FAILURE
    return
  endif

  ! Create src field
  call ESMF_ArraySpecSet(arrayspec1, 1, ESMF_TYPEKIND_R8, rc=rc)

  srcField = ESMF_FieldCreate(srcLocStream, arrayspec1, &
                              name="source", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif

!  if (petCount .eq. 1) then
    do lDE=0,localDECount-1
      call  ESMF_LocStreamGetBounds(srcLocStream, localDE=lDE, &
                              computationalLBound=cl, computationalUBound=cu, &
                              computationalCount=cc, &
                              rc=localrc)
      if (localrc /=ESMF_SUCCESS) then
        print*,'ERROR:  trouble with locStreamGet'
        rc=ESMF_FAILURE
        return
      endif


      !-------------------------------------------------------------------
      ! Get key data.
      !-------------------------------------------------------------------
      call ESMF_LocStreamGetKey(srcLocStream,                 &
                                localDE=lDE,                    &
                                keyName="ESMF:Lat",           &
                                farray=Yarray,                &
                                rc=localrc)
      if (localrc /=ESMF_SUCCESS) then
        print*,'ERROR:  trouble getting LocStream key for latitude'
        rc=ESMF_FAILURE
        return
      endif
      call ESMF_LocStreamGetKey(srcLocStream,                 &
                                localDE=lDE,                    &
                                keyName="ESMF:Lon",           &
                                farray=Xarray,                &
                                rc=localrc)
      if (localrc /=ESMF_SUCCESS) then
        print*,'ERROR:  trouble getting LocStream key for longitude'
        rc=ESMF_FAILURE
        return
      endif


      ! get src pointer
      call ESMF_FieldGet(srcField, lDE, srcPtr, rc=localrc)
      if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
      endif

      do idx=cl,cu
        i1=(idx-1)/src_nx + 1
        i2=mod((idx-1),src_nx) + 1

        ! Set source coordinates as 0 to 360
        Xarray(idx) = REAL(i2-1)*src_dx
        Yarray(idx) = -90. + (REAL(i1-1)*src_dy + 0.5*src_dy)

        ! Set the source to be a function of the x,y,z coordinate
        theta = DEG2RAD*(Xarray(idx))
        phi = DEG2RAD*(90.-Yarray(idx))

        srcPtr(idx) = (2. + cos(theta)**2.*cos(2.*phi))
      enddo
    enddo
!  endif



  ! setup dest. grid
  dstGrid=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


  call ESMF_ArraySpecSet(arrayspec2, 2, ESMF_TYPEKIND_R8, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif




  dstField = ESMF_FieldCreate(dstGrid, arrayspec2, &
                         staggerloc=ESMF_STAGGERLOC_CENTER, name="dest", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif


  xdstField = ESMF_FieldCreate(dstGrid, arrayspec2, &
                         staggerloc=ESMF_STAGGERLOC_CENTER, name="dest", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif



  call ESMF_GridAddCoord(dstGrid, staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif



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

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Destination grid
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! Get memory and set coords for dst
  do lDE=0,localDECount-1

     !! get coords
     call ESMF_GridGetCoord(dstGrid, 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(dstGrid, 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(dstField, lDE, dstPtr,  rc=localrc)
     if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
     endif


     call ESMF_FieldGet(xdstField, lDE, xdstPtr,  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)

        ! Set source coordinates as 0 to 360
        farrayPtrXC(i1,i2) = REAL(i1-1)*dst_dx
        farrayPtrYC(i1,i2) = -90. + (REAL(i2-1)*dst_dy + 0.5*dst_dy)

        ! initialize destination field
        dstPtr(i1,i2)=0.0

       ! Set the source to be a function of the x,y,z coordinate
        theta = DEG2RAD*(farrayPtrXC(i1,i2))
        phi = DEG2RAD*(90.-farrayPtrYC(i1,i2))

        ! After calculating field shift coords slighlty to be close, but not exact
        ! to make test more interesting
        farrayPtrXC(i1,i2) = farrayPtrXC(i1,i2) + 2.0

        xdstPtr(i1,i2) = (2. + cos(theta)**2.*cos(2.*phi))

     enddo
     enddo

  enddo    ! lDE


  ! Regrid store
  ! Calculate routeHandle on 2D fields
  call ESMF_FieldRegridStore( &
          srcField, &
          dstField=dstField, &
          routeHandle=routeHandle, &
          regridmethod=ESMF_REGRIDMETHOD_NEAREST_DTOS, &
          rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
      rc=ESMF_FAILURE
      return
   endif


  ! Do regrid on fields with extra dimension
  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


  ! Check results
  do lDE=0,localDECount-1

     ! Get interpolated dst field
     call ESMF_FieldGet(dstField, lDE, dstPtr, rc=localrc)
     if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
     endif

     call ESMF_FieldGet(xdstField, lDE, xdstPtr, computationalLBound=fclbnd, &
                             computationalUBound=fcubnd,  rc=localrc)
     if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
     endif

     ! Make sure everthing looks ok
     do i1=fclbnd(1),fcubnd(1)
     do i2=fclbnd(2),fcubnd(2)

        if (xdstPtr(i1,i2) .ne. 0.0) then
           if (abs(dstPtr(i1,i2)-xdstPtr(i1,i2))/abs(dstPtr(i1,i2)) &
                .gt. 0.05) then
              correct=.false.
              write(*,*) i1,i2,"::",dstPtr(i1,i2),xdstPtr(i1,i2),(dstPtr(i1,i2)-xdstPtr(i1,i2))/xdstPtr(i1,i2)
           endif
        else
           if (abs(dstPtr(i1,i2)-xdstPtr(i1,i2)) &
                .gt. 0.05) then
              correct=.false.
           endif
        endif
     enddo
     enddo

  enddo    ! lDE


  ! 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

   call ESMF_FieldDestroy(xdstField, rc=localrc)
   if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
   endif


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

  call ESMF_LocStreamDestroy(srcLocStream, 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_regridNearestLocStreamToGrid