test_regridGridToLocStreamRegDist Subroutine

subroutine test_regridGridToLocStreamRegDist(rc)

Arguments

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

Source Code

 subroutine test_regridGridToLocStreamRegDist(rc)
  integer, intent(out)  :: rc
  logical :: correct
  integer :: localrc
  type(ESMF_Grid) :: srcGrid
  type(ESMF_LocStream) :: dstLocStream
  type(ESMF_Field) :: srcField,dstField,dstFieldPatch
  type(ESMF_RouteHandle) :: routeHandle,routeHandlePatch
  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(:,:), farrayPtr1D(:)
  real(ESMF_KIND_R8), pointer :: latArray(:),lonArray(:)
  real(ESMF_KIND_R8) :: src_dx, src_dy, dst_dx, dst_dy
  real(ESMF_KIND_R8) :: RAD2DEG,DEG2RAD,theta,phi
  real(ESMF_KIND_R8) :: lat,lon
  real(ESMF_KIND_R8) :: x,y,z,expected
  integer :: clbnd(2),cubnd(2)
  integer :: fclbnd(2),fcubnd(2)
  integer :: dclbnd(1),dcubnd(1)
  integer :: i1,i2
  integer :: lDE, localDECount
  integer :: src_nx, src_ny, dst_nx, dst_ny
  integer :: cl,cu,idx
  integer :: localPet, petCount
  real(ESMF_KIND_R8) :: beg_time, end_time
#if defined (ESMF_LAPACK)
  logical, external :: LSAME
  logical :: tf
#endif

  ! 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

  ! If we don't have 1 or 4 PETS then exit with failure
  if ((petCount .ne. 1) .and. (petCount .ne. 4)) then
    print*,'ERROR:  test must be run using exactly 1 or 4 PETS - detected ',petCount
    rc=ESMF_FAILURE
    return
   endif


  ! Establish the resolution of the grids
  src_nx = 100
  src_ny = 50
  src_dx = 360./src_nx
  src_dy = 180./src_ny
  DEG2RAD = 3.14159265/180.0
  RAD2DEG = 1./DEG2RAD

  ! setup source grid
  srcGrid=ESMF_GridCreate1PeriDim(minIndex=(/1,1/),maxIndex=(/src_nx,src_ny/),regDecomp=(/petCount,1/), &
                                  coordSys=ESMF_COORDSYS_SPH_DEG, indexflag=ESMF_INDEX_GLOBAL, &
                                  rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif

  ! Create source field
  call ESMF_ArraySpecSet(arrayspec, 2, ESMF_TYPEKIND_R8, rc=rc)

  srcField = ESMF_FieldCreate(srcGrid, arrayspec, &
                              staggerloc=ESMF_STAGGERLOC_CENTER, name="source", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif

 ! Allocate coordinates
  call ESMF_GridAddCoord(srcGrid, staggerloc=ESMF_STAGGERLOC_CENTER, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    rc=ESMF_FAILURE
    return
  endif

  ! 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 src
  do lDE=0,localDECount-1

     !! get coord 1
     call ESMF_GridGetCoord(srcGrid, 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(srcGrid, 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(srcField, 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 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

     enddo
     enddo

  enddo    ! lDE




  ! Setup Dst LocStream
  dst_nx = 90
  dst_ny = 40
  dst_dx = 360./dst_nx
  dst_dy = 180./dst_ny

  dstLocStream=ESMF_LocStreamCreate(minIndex=1, maxIndex=dst_nx*dst_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(dstLocStream, 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(dstLocStream,                 &
                            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(dstLocStream,                 &
                            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


  !-------------------------------------------------------------------
  ! Get key data.
  !-------------------------------------------------------------------
  do lDE=0,localDECount-1
    call  ESMF_LocStreamGetBounds(dstLocStream, localDE=lDE, &
                            computationalLBound=cl, computationalUBound=cu, &
                            rc=localrc)
    if (localrc /=ESMF_SUCCESS) then
      print*,'ERROR:  trouble with LocStreamGet'
      rc=ESMF_FAILURE
      return
    endif

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

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

      ! Set source coordinates as 0 to 360
      lonArray(idx) = REAL(i2-1)*dst_dx
      latArray(idx) = -90. + (REAL(i1-1)*dst_dy + 0.5*dst_dy)
    enddo
  enddo

  ! Create dest fields
  call ESMF_ArraySpecSet(arrayspec, 1, ESMF_TYPEKIND_R8, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    print*,'ERROR:  trouble calling ArraySpecSet'
    rc=ESMF_FAILURE
    return
   endif

  dstField = ESMF_FieldCreate(dstLocStream, arrayspec, &
                        name="dest", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    print*,'ERROR:  trouble creating field on locStream'
    rc=ESMF_FAILURE
    return
  endif

  dstFieldPatch = ESMF_FieldCreate(dstLocStream, arrayspec, &
                        name="destPatch", rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
    print*,'ERROR:  trouble creating field on locStream'
    rc=ESMF_FAILURE
    return
  endif

  ! clear destination Fields
  ! Should only be 1 localDE
  call ESMF_FieldGet(dstField, 0, farrayPtr1D,  rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
  endif

  farrayPtr1D=0.0

  call ESMF_FieldGet(dstFieldPatch, 0, farrayPtr1D,  rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
  endif

  farrayPtr1D=0.0

  !!! Regrid forward from the grid to the LocStream
  ! 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

  !!! Regrid forward from the Src grid to the LocStream - this time with PATCH

  ! Make sure LAPACK lib is mapped in before the timing block
#if defined (ESMF_LAPACK)
  tf = LSAME ('a', 'A')
#endif

  ! Get start time
  call ESMF_VMWtime(beg_time)

  ! Regrid store
  call ESMF_FieldRegridStore( &
          srcField, &
          dstField=dstFieldPatch, &
          routeHandle=routeHandlePatch, &
          regridmethod=ESMF_REGRIDMETHOD_PATCH, &
          rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
      rc=localrc
      return
   endif

  ! Get end time
  call ESMF_VMWtime(end_time)
  write(*,*) localPet," Time to do Patch FieldRegridStore()=",end_time-beg_time

  ! Do regrid
  call ESMF_FieldRegrid(srcField, dstFieldPatch, routeHandlePatch, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
      rc=ESMF_FAILURE
      return
   endif

   call ESMF_FieldRegridRelease(routeHandlePatch, rc=localrc)
  if (localrc /=ESMF_SUCCESS) then
      rc=ESMF_FAILURE
      return
   endif

  ! Check results
  do lDE=0,localDECount-1

     call ESMF_FieldGet(dstField, lDE, farrayPtr1D, computationalLBound=dclbnd, &
                             computationalUBound=dcubnd,  rc=localrc)
     if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
     endif

     ! Make sure everthing looks ok
     do i1=dclbnd(1),dcubnd(1)
        ! Get coordinates
        lon=lonArray(i1)
        lat=latArray(i1)

        ! get the x,y,z coordinates
        theta = DEG2RAD*(lon)
        phi = DEG2RAD*(90.-lat)
        x = cos(theta)*sin(phi)
        y = sin(theta)*sin(phi)
        z = cos(phi)

        ! determine validation data
        expected = x+y+z+15.0

        !! if error is too big report an error
        if ( abs( farrayPtr1D(i1)-(expected) )/expected > 0.001) then
           print*,'ERROR: larger than expected difference, expected ',expected, &
                  '  got ',farrayPtr1D(i1),'  diff= ',abs(farrayPtr1D(i1)-expected)
           correct=.false.
        endif
     enddo

     ! now for patch
     call ESMF_FieldGet(dstFieldPatch, lDE, farrayPtr1D, computationalLBound=dclbnd, &
                             computationalUBound=dcubnd,  rc=localrc)
     if (localrc /=ESMF_SUCCESS) then
        rc=ESMF_FAILURE
        return
     endif

     ! Make sure everthing looks ok
     do i1=dclbnd(1),dcubnd(1)
        ! Get coordinates
        lon=lonArray(i1)
        lat=latArray(i1)

        ! get the x,y,z coordinates
        theta = DEG2RAD*(lon)
        phi = DEG2RAD*(90.-lat)
        x = cos(theta)*sin(phi)
        y = sin(theta)*sin(phi)
        z = cos(phi)

        ! determine validation data
        expected = x+y+z+15.0

        !! if error is too big report an error
        if ( abs( farrayPtr1D(i1)-(expected) )/expected > 0.001) then
           print*,'ERROR: larger than expected difference, expected ',expected, &
                  '  got ',farrayPtr1D(i1),'  diff= ',abs(farrayPtr1D(i1)-expected)
           correct=.false.
        endif
     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(dstFieldPatch, rc=localrc)
   if (localrc /=ESMF_SUCCESS) then
     rc=ESMF_FAILURE
     return
   endif

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

  call ESMF_GridDestroy(srcGrid, 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_regridGridToLocStreamRegDist