computeAreaGrid Subroutine

public subroutine computeAreaGrid(grid, petNo, area, regridScheme, rc)

Arguments

Type IntentOptional Attributes Name
type(ESMF_Grid) :: grid
integer :: petNo
real(kind=ESMF_KIND_R8), pointer :: area(:)
integer :: regridScheme
integer :: rc

Source Code

subroutine computeAreaGrid(grid, petNo, area, regridScheme, rc)
  type(ESMF_Grid) :: grid
  integer :: petNo
  real (ESMF_KIND_R8), pointer :: area(:)
  integer :: regridScheme
  integer :: rc

  type(ESMF_Field) :: areaField
  integer :: minIndex(2), maxIndex(2), gridDims(2)
  real (ESMF_KIND_R8), pointer :: area2D(:,:)
  integer :: i, start, ntiles
  integer :: localrc

  ! Create a field on the grid to hold the areas
  areaField=ESMF_FieldCreate(grid, typekind=ESMF_TYPEKIND_R8, &
                 staggerloc=ESMF_STAGGERLOC_CENTER, name="area", rc=localrc)
  if (ESMF_LogFoundError(localrc, &
                         ESMF_ERR_PASSTHRU, &
                         ESMF_CONTEXT, rcToReturn=rc)) return
  if (localrc /=ESMF_SUCCESS) then
     rc=localrc
     return
  endif

  ! compute areas
  call ESMF_FieldRegridGetArea(areaField, rc=localrc)
  if (ESMF_LogFoundError(localrc, &
                         ESMF_ERR_PASSTHRU, &
                         ESMF_CONTEXT, rcToReturn=rc)) return
 if (localrc /=ESMF_SUCCESS) then
     rc=localrc
     return
  endif

  ! Get number of tiles
  call ESMF_GridGet(grid, tileCount = ntiles, rc=localrc)
  if (ESMF_LogFoundError(localrc, &
       ESMF_ERR_PASSTHRU, &
       ESMF_CONTEXT, rcToReturn=rc)) return

  ! Get size of Grid
  call ESMF_GridGet(grid, tile=1, staggerloc=ESMF_STAGGERLOC_CENTER, &
        minIndex=minIndex, maxIndex=maxIndex, rc=localrc)
  if (ESMF_LogFoundError(localrc, &
                         ESMF_ERR_PASSTHRU, &
                         ESMF_CONTEXT, rcToReturn=rc)) return
  if (localrc /=ESMF_SUCCESS) then
      rc=localrc
      return
  endif

  ! Grid size
  gridDims(1)=maxIndex(1)-minIndex(1)+1
  gridDims(2)=maxIndex(2)-minIndex(2)+1

  ! Allocate memory for area
  allocate(area2D(gridDims(1),gridDims(2)))

  ! Only do this part on PET 0
  if (petNo .eq. 0) then
     ! Allocate memory for area
     allocate(area(gridDims(1)*gridDims(2)*ntiles))
  endif

  ! Get area onto PET 0
  start=1
  do i=1,ntiles
    call ESMF_FieldGather(areaField, farray=area2D, rootPet=0, tile=i, rc=localrc)
    if (ESMF_LogFoundError(localrc, &
                         ESMF_ERR_PASSTHRU, &
                         ESMF_CONTEXT, rcToReturn=rc)) return

     ! copy to 1D array
     if (PetNo == 0) then
       ! flatten area
       area(start:start+gridDims(1)*gridDims(2)-1)=RESHAPE(area2D,(/gridDims(1)*gridDims(2)/))
       start= start+gridDims(1)*gridDims(2)
     endif
  enddo

  ! deallocate memory for 2D area
  deallocate(area2D)

end subroutine computeAreaGrid