areamod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

! Get indices and weights for mapping from input grid to output grid
! -----------------------------------------------------------------

    call areamap_point (io      , jo       , nlon_i  , nlat_i , nlon_o  , &
                        nlat_o  , numlon_i , lon_i   , lat_i  , mask_i  , &
                        lon_o   , lat_o    , fland_o , area_o , novr_i2o, &
                        iovr_i2o, jovr_i2o , wovr_i2o, lon_i_offset)

! -----------------------------------------------------------------
! Error check domain 
! -----------------------------------------------------------------

    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
       write (6,*) 'AREAINI warning: conservation check 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

    return
  end subroutine areaini_point

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

  subroutine areamap_point (io     , jo       , nlon_i  , nlat_i , nlon_o , &
                            nlat_o , numlon_i , lon_i   , lat_i  , mask_i , &
                            lon_o  , lat_o    , fland_o , area_o , n_ovr  , &
                            i_ovr  , j_ovr    , w_ovr   , lon_i_offset)

!----------------------------------------------------------------------- 
! 
! 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,maxovr),j_ovr(io,jo,maxovr)) * w_ovr(io,jo,maxovr) 
!
! Note: maxovr is some number greater than n_ovr. Weights of zero are 
! used for the excess points
! 
! Author: Gordon Bonan
! 
!-----------------------------------------------------------------------

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

! ------------------------ local variables ---------------------------
    integer  :: ii                  !input  grid longitude 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
! --------------------------------------------------------------------

    n_ovr           = 0
    i_ovr(1:maxovr) = 1
    j_ovr(1:maxovr) = 1
    w_ovr(1:maxovr) = 0. 

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

    call areaovr_point (io    , jo    , nlon_i , nlat_i , numlon_i , &
                        lon_i , lat_i , nlon_o , nlat_o , lon_o    , &
                        lat_o , n_ovr , i_ovr  , j_ovr  , w_ovr    )

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

    call areaovr_point (io           , jo    , nlon_i , nlat_i , numlon_i , &
                        lon_i_offset , lat_i , nlon_o , nlat_o , lon_o    , &
                        lat_o        , n_ovr , i_ovr  , j_ovr  , w_ovr    )

! --------------------------------------------------------------------
! 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.
! --------------------------------------------------------------------

! find total land area of overlapping input cells

    f_ovr = 0.
    do n = 1, n_ovr
       ii = i_ovr(n)
       ji = j_ovr(n)
       f_ovr = f_ovr + w_ovr(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)/area_o > 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
       call endrun
    end if

! make weights

    do n = 1, n_ovr
       ii = i_ovr(n)
       ji = j_ovr(n)
       if (f_ovr > 0.) then
          w_ovr(n) = w_ovr(n)*mask_i(ii,ji) / f_ovr
       else
          w_ovr(n) = 0.
       end if
    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

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

    f_ovr = 0.
    do n = 1, maxovr
       f_ovr = f_ovr + w_ovr(n)
    end do

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

    return
  end subroutine areamap_point

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

  subroutine areaovr_point (io     , jo     , nlon_i , nlat_i , numlon_i , &
                            lon_i  , lat_i  , nlon_o , nlat_o , lon_o    , &
                            lat_o  , 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
!
!           lon_o(io,jo)      lon_o(io+1,jo)
!
!              |                   |
!              --------------------- lat_o(jo+1)
!              |                   |
!              |                   |
!    xxxxxxxxxxxxxxx lat_i(ji+1)   |
!    x         |   x               |
!    x  input  |   x   output      |
!    x  cell   |   x    cell       |
!    x  ii,ji  |   x   io,jo       |
!    x         |   x               |
!    x         ----x---------------- lat_o(jo  )
!    x             x
!    xxxxxxxxxxxxxxx lat_i(ji  )
!    x             x
! lon_i(ii,ji) lon_i(ii+1,ji)
!
!
! The above diagram assumes both grids are oriented South to North. Other
! combinations of North to South and South to North grids are possible:
!
!     Input Grid    Output Grid
!     -------------------------
! (1)   S to N        S to N
! (2)   N to S        N to S
! (3)   S to N        N to S
! (4)   N to S        S to N
!
! The code has been modified to allow for North to South grids. Verification
! that these changes work are: 
!    o (1) and (4) give same results for output grid
!    o (2) and (3) give same results for output grid
!    o (2) and (4) give same results for output grid when output grid inverted
!
! WARNING: this code does not vectorize but is only called during start-up
! 
! Author: Gordon Bonan
! 
!-----------------------------------------------------------------------

! ------------------------ arguments -----------------------------------
    integer , intent(in) :: io                     !output grid lon index 
    integer , intent(in) :: jo                     !output grid lat index   
    integer , intent(in) :: nlon_i                 !input grid : max number of longitude points
    integer , int