areamod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

       end do
    end do
!$OMP END PARALLEL DO

! loop through overlapping cells on input grid to make area-average

    do n = 1, nmax
!$OMP PARALLEL DO PRIVATE (jo,io,ii,ji)
       do jo = 1, nlat_o
          do io =1, numlon_o(jo)
             ii = i_ovr(io,jo,n)
             ji = j_ovr(io,jo,n)
             fld_o(io,jo) = fld_o(io,jo) + w_ovr(io,jo,n)*fld_i(ii,ji)
          end do
       end do
!$OMP END PARALLEL DO
    end do

! set non-valid points for reduced grid to missing value

    do jo = 1, nlat_o
       do io = numlon_o(jo)+1, nlon_o
          fld_o(io,jo) = 1.e36
       end do
    end do

    return
  end subroutine areaave

!=======================================================================

  subroutine areamap (nlon_i   , nlat_i   , nlon_o , nlat_o ,  &
                      lon_i    , lat_i    , lon_o  , lat_o  ,  &
                      numlon_i , numlon_o , mask_i , mx_ovr ,  &
                      n_ovr    , i_ovr    , j_ovr  , w_ovr  ,  &
                      fland_o  , area_o   )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! weights and indices for area of overlap between grids
!
! Method: 
! Get indices and weights for area-averaging between input and output grids.
! For each output grid cell find:
!
!    o number of input grid cells that overlap with output grid cell (n_ovr)
!    o longitude index (1 <= i_ovr <= nlon_i) of the overlapping input grid cell
!    o latitude index  (1 <= j_ovr <= nlat_i) of the overlapping input grid cell
!    o fractional overlap of input grid cell (w_ovr)
!
! so that for
!
! field values fld_i on an  input grid with dimensions nlon_i and nlat_i
! field values fld_o on an output grid with dimensions nlon_o and nlat_o are
!
! fld_o(io,jo) = 
! fld_i(i_ovr(io,jo,     1),j_ovr(io,jo,     1)) * w_ovr(io,jo,     1) +
!                             ... + ... +
! fld_i(i_ovr(io,jo,mx_ovr),j_ovr(io,jo,mx_ovr)) * w_ovr(io,jo,mx_ovr) 
!
! Note: mx_ovr is some number greater than n_ovr. Weights of zero are 
! used for the excess points
! 
! Author: Gordon Bonan
! 
!-----------------------------------------------------------------------

! ------------------------ arguments ---------------------------------
    integer ,intent(in) :: nlon_i                     !input grid : max number of longitude points
    integer ,intent(in) :: nlat_i                     !input grid : number of latitude points
    integer ,intent(in) :: numlon_i(nlat_i)           !input grid : number longitude points for lat
    real(r8),intent(inout) :: lon_i(nlon_i+1,nlat_i)  !input grid : cell longitude, west edge (deg)
    real(r8),intent(in) :: lat_i(nlat_i+1)            !input grid : cell latitude,  south edge (deg)
    real(r8),intent(in) :: mask_i(nlon_i,nlat_i)      !input grid : mask (0, 1) 
    integer ,intent(in) :: nlon_o                     !output grid: max number of longitude points
    integer ,intent(in) :: nlat_o                     !output grid: number of latitude points
    integer ,intent(in) :: numlon_o(nlat_o)           !output grid: number longitude points for lat
    real(r8),intent(in) :: lon_o(nlon_o+1,nlat_o)     !output grid: cell longitude, west edge  (deg)
    real(r8),intent(in) :: lat_o(nlat_o+1)            !output grid: cell latitude,  south edge (deg)
    real(r8),intent(in) :: fland_o(nlon_o,nlat_o)     !output grid: fraction that is land
    real(r8),intent(in) :: area_o(nlon_o,nlat_o)      !output grid: cell area
    integer ,intent(in) :: mx_ovr                     !max num input cells that overlap output cell
    integer ,intent(out):: n_ovr(nlon_o,nlat_o)       !number of overlapping input cells
    integer ,intent(out):: i_ovr(nlon_o,nlat_o,mx_ovr)!lon index, overlapping input cell
    integer ,intent(out):: j_ovr(nlon_o,nlat_o,mx_ovr)!lat index, overlapping input cell
    real(r8),intent(out):: w_ovr(nlon_o,nlat_o,mx_ovr)!overlap weights for input cells
! --------------------------------------------------------------------

! ------------------------ local variables ---------------------------
    integer :: io                   !output grid longitude loop index
    integer :: ii                   !input  grid longitude loop index
    integer :: jo                   !output grid latitude  loop index
    integer :: ji                   !input  grid latitude  loop index
    integer :: n                    !overlapping cell index
    real(r8) :: offset              !used to shift x-grid 360 degrees
    real(r8) :: f_ovr               !sum of overlap weights 
    real(r8) :: relerr = 0.00001    !max error: sum overlap weights ne 1
    real(r8) :: dx_i                !input grid  longitudinal range
    real(r8) :: dy_i                !input grid  latitudinal  range
    real(r8) :: dx_o                !output grid longitudinal range
    real(r8) :: dy_o                !output grid latitudinal  range
! --------------------------------------------------------------------

! --------------------------------------------------------------------
! Initialize overlap weights on output grid to zero for maximum 
! number of overlapping points. Set lat and lon indices of overlapping
! input cells to dummy values. Set number of overlapping cells to zero
! --------------------------------------------------------------------

    do n = 1, mx_ovr 
       do jo = 1, nlat_o      
          do io = 1, numlon_o(jo)  
             i_ovr(io,jo,n) = 1
             j_ovr(io,jo,n) = 1
             w_ovr(io,jo,n) = 0. 
          end do
       end do
    end do

    do jo = 1, nlat_o      
       do io = 1, numlon_o(jo)  
          n_ovr(io,jo) = 0
       end do
    end do

! --------------------------------------------------------------------
! First pass to find cells that overlap and area of overlap
! --------------------------------------------------------------------

    call areaovr (nlon_i , nlat_i , numlon_i , lon_i  , lat_i  , &
                  nlon_o , nlat_o , numlon_o , lon_o  , lat_o  , &
                  mx_ovr , n_ovr  , i_ovr    , j_ovr  , w_ovr  )

! --------------------------------------------------------------------
! Second pass to find cells that overlap and area of overlap
! --------------------------------------------------------------------

! Shift x-grid to locate periodic grid intersections. This
! assumes that all lon_i(1,j) have the same value for all
! latitudes j and that the same holds for lon_o(1,j)

    if (lon_i(1,1) < lon_o(1,1)) then
       offset = 360.0
    else
       offset = -360.0
    end if

    do ji = 1, nlat_i
       do ii = 1, numlon_i(ji) + 1
          lon_i(ii,ji) = lon_i(ii,ji) + offset
       end do
    end do

! find overlap

    call areaovr (nlon_i , nlat_i , numlon_i , lon_i  , lat_i  , &
                  nlon_o , nlat_o , numlon_o , lon_o  , lat_o  , &
                  mx_ovr , n_ovr  , i_ovr    , j_ovr  , w_ovr  )

! restore x-grid (un-shift x-grid)

    do ji = 1, nlat_i
       do ii = 1, numlon_i(ji) + 1
          lon_i(ii,ji) = lon_i(ii,ji) - offset
       end do
    end do

! --------------------------------------------------------------------
! Normalize areas of overlap to get fractional contribution of each
! overlapping grid cell (input grid) to grid cell average on output grid.
! Normally, do this by dividing area of overlap by area of output grid cell.
! But, only have data for land cells on input grid. So if output grid cell
! overlaps with land and non-land cells (input grid), do not have valid
! non-land data for area-average. Instead, weight by area of land using
! [mask_i], which has a value of one for land and zero for ocean. If
! [mask_i] = 1, input grid cell contributes to output grid cell average.
! If [mask_i] = 0, input grid cell does not contribute to output grid cell 
! average.
! --------------------------------------------------------------------

    do jo = 1, nlat_o
       do io = 1, numlon_o(jo)

! find total land area of overlapping input cells

          f_ovr = 0.
          do n = 1, n_ovr(io,jo)
             ii = i_ovr(io,jo,n)
             ji = j_ovr(io,jo,n)
             f_ovr = f_ovr + w_ovr(io,jo,n)*mask_i(ii,ji)
          end do

! make sure area of overlap is less than or equal to output grid cell area

          if ((f_ovr-area_o(io,jo))/area_o(io,jo) > relerr) then
             write (6,*) 'AREAMAP error: area not conserved for lon,lat = ',io,jo
             write (6,'(a30,e20.10)') 'sum of overlap area = ',f_ovr
             write (6,'(a30,e20.10)') 'area of output grid = ',area_o(io,jo)
             call endrun
          end if

! make weights

          do n = 1, n_ovr(io,jo)
             ii = i_ovr(io,jo,n)
             ji = j_ovr(io,jo,n)
             if (f_ovr > 0.) then
                w_ovr(io,jo,n) = w_ovr(io,jo,n)*mask_i(ii,ji) / f_ovr
             else
                w_ovr(io,jo,n) = 0.
             end if
          end do

       end do
    end do

! --------------------------------------------------------------------
! Error check: overlap weights for input grid cells must sum to 1. This
! is always true if both grids span the same domain. However, if one
! grid is a subset of the other grid, this is only true when mapping
! from the full grid to the subset. When input grid covers a smaller
! domain than the output grid, this test is not valid.
! --------------------------------------------------------------------

    dx_i = lon_i(nlon_i+1,1) - lon_i(1,1)
    dx_o = lon_o(nlon_o+1,1) - lon_o(1,1)

    if (lat_i(nlat_i+1) > lat_i(1)) then      !South to North grid
       dy_i = lat_i(nlat_i+1) - lat_i(1)
    else                                         !North to South grid
       dy_i = lat_i(1) - lat_i(nlat_i+1)
    end if
    if (lat_o(nlat_o+1) > lat_o(1)) then      !South to North grid
       dy_o = lat_o(nlat_o+1) - lat_o(1)
    else                                         !North to South grid
       dy_o = lat_o(1) - lat_o(nlat_o+1)
    end if

    if (abs(dx_i-dx_o)>relerr .or. abs(dy_i-dy_o)>relerr) then
       if (dx_i<dx_o .or. dy_i<dy_o) then
          write (6,*) 'AREAMAP warning: area-average not valid for '
          write (6,*) '   input  grid of ',nlon_i,' x ',nlat_i
          write (6,*) '   output grid of ',nlon_o,' x ',nlat_o
          return     
       end if
    end if

    do jo = 1, nlat_o
       do io = 1, numlon_o(jo)
          f_ovr = 0.

          do n = 1, mx_ovr
             f_ovr = f_ovr + w_ovr(io,jo,n)
          end do

! error check only valid if output grid cell has land. non-land cells
! will have weights equal to zero

          if (fland_o(io,jo) > 0.) then
             if (abs(f_ovr-1.) > relerr) then
                write (6,*) 'AREAMAP error: area not conserved for lon,lat = ',io,jo
                write (6,'(a30,e20.10)') 'sum of overlap weights = ',f_ovr
                call endrun
             end if
          end if

       end do
    end do

    return
  end subroutine areamap

!=======================================================================

  subroutine areaovr (nlon_i , nlat_i , numlon_i , lon_i , lat_i  , &
                      nlon_o , nlat_o , numlon_o , lon_o , lat_o  , &
                      mx_ovr , n_ovr  , i_ovr    , j_ovr , w_ovr  )


!----------------------------------------------------------------------- 
! 
! Purpose: 
! Find area of overlap between grid cells
!
! Method: 
! For each output grid cell: find overlapping input grid cell and area of
! input grid cell that overlaps with output grid cell. Cells overlap if:
!
! southern edge of input grid < northern edge of output grid AND
! northern edge of input grid > southern edge of output grid
!
! western edge of input grid < eastern edge of output grid AND
! eastern edge of input grid > western edge of output grid
!