areamod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

! Method: 
! 
! Author: Mariana Vertenstein
! 
!-----------------------------------------------------------------------

! ------------------------ arguments ---------------------------------
    integer , intent(in) :: nlat              !dimension: number of latitude points
    integer , intent(in) :: nlon              !dimension: number of longitude points
    integer , intent(in) :: numlon(nlat)      !number of grid cells per latitude strip
    real(r8), intent(in) :: longxy(nlon,nlat) !longitude at center of grid cell
    real(r8), intent(in) :: latixy(nlon,nlat) !latitude at center of grid cell
    real(r8), intent(in) :: edgen             !northern edge of grid (degrees)
    real(r8), intent(in) :: edgee             !eastern edge of grid (degrees)
    real(r8), intent(in) :: edges             !southern edge of grid (degrees)
    real(r8), intent(in) :: edgew             !western edge of grid (degrees)
    real(r8), intent(out):: lats(nlat+1)      !grid cell latitude, southern edge (degrees)
    real(r8), intent(out):: lonw(nlon+1,nlat) !grid cell longitude, western edge (degrees)
! --------------------------------------------------------------------

! ------------------- local variables -----------------------------
    integer i,j             !indices
    real(r8) dx             !cell width
! --------------------------------------------------------------------

! --------------------------------------------------------------------
! Latitudes -- southern/northern edges for each latitude strip. 
! For grids oriented South to North, the southern
! and northern edges of latitude strip [j] are:
!        southern = lats(j  )
!        northern = lats(j+1)
! For grids oriented North to South: the southern
! and northern edges of latitude strip [j] are:
!        northern = lats(j  )
!        southern = lats(j+1)
! In both cases, [lats] must be dimensioned lats(lat+1)
! --------------------------------------------------------------------

    if (nlat == 1) then
       lats(1) = edges 
       lats(nlat+1) = edgen 
    else
       if (latixy(1,2) > latixy(1,1)) then    !South to North grid
          lats(1) = edges 
          lats(nlat+1) = edgen 
       else                                   !North to South grid
          lats(1) = edgen 
          lats(nlat+1) = edges 
       end if
    end if
    do j = 2, nlat
       lats(j) = (latixy(1,j-1) + latixy(1,j)) / 2.
    end do

! --------------------------------------------------------------------
! Longitudes -- western edges. Longitudes for the western edge of the 
! cells must increase continuously and span 360 degrees. Assume that
! grid starts at Dateline with western edge on Dateline Western edges 
! correspond to [longxy] (longitude at center of cell) and range from 
! -180 to 180 with negative longitudes west of Greenwich. 
! Partial grids that do not span 360 degrees are allowed so long as they
! have the convention of Grid 1 with 
!      western edge of grid: >= -180 and < 180
!      eastern edge of grid: > western edge  and <= 180
! [lonw] must be dimensioned lonw(lon+1,lat) because each latitude
! strip can have variable longitudinal resolution
! --------------------------------------------------------------------

! Western edge of first grid cell -- since grid starts with western
! edge on Dateline, lonw(1,j)=-180. This is the same as [edgew].
! Remaining grid cells. On a global grid lonw(numlon+1,j)=lonw(1,j)+360. 
! This is the same as [edgee].  Set unused longitudes to non-valid number

    do j = 1, nlat
       dx = (edgee - edgew) / numlon(j)
       lonw(1,j) = edgew
       do i = 2, numlon(j)+1
          lonw(i,j) = lonw(1,j) + (i-1)*dx
       end do
       do i = numlon(j)+2, nlon
          lonw(i,j) = -999.
       end do
    end do

    return
  end subroutine celledge_regional

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

  subroutine celledge_global (nlat, nlon, numlon, longxy, latixy, &
                              lats, lonw )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Southern and western edges of grid cells - global grid
!
! Method: 
! 
! Author: Mariana Vertenstein
! 
!-----------------------------------------------------------------------

! ------------------------ arguments ---------------------------------
    integer , intent(in) :: nlat               !dimension: number of latitude points
    integer , intent(in) :: nlon               !dimension: number of longitude points
    integer , intent(in) :: numlon(nlat)       !number of grid cells per latitude strip
    real(r8), intent(in) :: longxy(nlon,nlat)  !longitude at center of grid cell
    real(r8), intent(in) :: latixy(nlon,nlat)  !latitude at center of grid cell
    real(r8), intent(out):: lats(nlat+1)       !grid cell latitude, southern edge (degrees)
    real(r8), intent(out):: lonw(nlon+1,nlat)  !grid cell longitude, western edge (degrees)
! --------------------------------------------------------------------

! ------------------- local variables -----------------------------
    integer i,j             !indices
    real(r8) dx             !cell width
! --------------------------------------------------------------------

! --------------------------------------------------------------------
! Latitudes -- southern/northern edges for each latitude strip. 
! For grids oriented South to North, the southern
! and northern edges of latitude strip [j] are:
!        southern = lats(j  )
!        northern = lats(j+1)
! For grids oriented North to South: the southern
! and northern edges of latitude strip [j] are:
!        northern = lats(j  )
!        southern = lats(j+1)
! In both cases, [lats] must be dimensioned lats(lat+1)
! --------------------------------------------------------------------

    do j = 1, nlat+1              !southern edges
       if (j == 1) then           !south pole
          lats(j) = -90.
       else if (j == nlat+1) then !north pole
          lats(j) = 90.
       else                       !edge = average latitude
          lats(j) = (latixy(1,j-1) + latixy(1,j)) / 2.
       end if
    end do

! --------------------------------------------------------------------
! Longitudes -- western edges. Longitudes for the western edge of the 
! cells must increase continuously and span 360 degrees.
! --------------------------------------------------------------------

    if (longxy(1,1) >= 0.) then
       do j = 1, nlat
          dx = 360./(numlon(j))
          do i = 1, numlon(j)+1
             lonw(i,j) = -dx/2. + (i-1)*dx
          end do
          do i = numlon(j)+2, nlon
             lonw(i,j) = -999.
          end do
       end do
    else
       write(6,*)'global non-regional grids currently only supported ', &
            'for grids starting at greenwich and centered on Greenwich'
       call endrun
    endif

    return
  end subroutine celledge_global

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

  subroutine areaini_point (io      , jo          , nlon_i  , nlat_i  , numlon_i, &
                            lon_i   , lon_i_offset, lat_i   , area_i  , mask_i  , &
                            nlon_o  , nlat_o      , numlon_o, lon_o   , lat_o   , &
                            area_o  , fland_o     , novr_i2o, iovr_i2o, jovr_i2o, &
                            wovr_i2o)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! area averaging initialization 
!
! Method: 
! This subroutine is used in conjunction with areaave.F for area-average
! mapping of a field from one grid to another. 
!
!    areaini_point  - initializes indices and weights for area-averaging from 
!                     input grid to output grid
!    areamap_point  - called by areaini: finds indices and weights
!    areaovr_onit   - called by areamap: finds if cells overlap and area of overlap
! 
! To map from one grid to another, must first call areaini to build
! the indices and weights (iovr_i2o, jovr_i2o, wovr_i2o). Then must
! call areaave to get new field on output grid.
!
! Not all grid cells on the input grid will be used in the area-averaging
! of a field to the output grid. Only input grid cells with [mask_i] = 1
! contribute to output grid cell average. If [mask_i] = 0, input grid cell 
! does not contribute to output grid cell. This distinction is not usually
! required for atm -> land mapping, because all cells on the atm grid have
! data. But when going from land -> atm, only land grid cells have data.
! Non-land grid cells on surface grid do not have data. So if output grid cell
! overlaps with land and non-land cells (input grid), can only use land
! grid cells when computing area-average. 
!
! o Input and output grids can be ANY resolution BUT:
!
!   a. Grid orientation -- Grids can be oriented south to north
!      (i.e. cell(lat+1) is north of cell(lat)) or from north to 
!      south (i.e. cell(lat+1) is south of cell(lat)). Both grids must be 
!      oriented from west to east, i.e., cell(lon+1) must be east of cell(lon)
!
!   b. Grid domain -- Grids do not have to be global. Both grids are defined
!      by their north, east, south, and west edges (edge_i and edge_o in
!      this order, i.e., edge_i(1) is north and edge_i(4) is west).
!      
!      For partial grids, northern and southern edges are any latitude
!      between 90 (North Pole) and -90 (South Pole). Western and eastern
!      edges are any longitude between -180 and 180, with longitudes 
!      west of Greenwich negative.
!
!      For global grids, northern and southern edges are 90 (North Pole) 
!      and -90 (South Pole). The grids do not have to start at the
!      same longitude, i.e., one grid can start at Dateline and go east;
!      the other grid can start at Greenwich and go east. Longitudes for
!      the western edge of the cells must increase continuously and span
!      360 degrees. Examples
!
!                              West edge    East edge
!                            ---------------------------------------------------
!      Dateline            :        -180 to 180        (negative W of Greenwich)
!      Greenwich (centered):    0 - dx/2 to 360 - dx/2 
!
!   c. Both grids can have variable number of longitude points for each
!      latitude strip. However, the western edge of the first point in each
!      latitude must be the same for all latitudes. Likewise, for the
!      eastern edge of the last point. That is, each latitude strip must span 
!      the same longitudes, but the number of points to do this can be different
!
!   d. One grid can be a sub-set (i.e., smaller domain) than the other grid.
!      In this way, an atmospheric dataset for the entire globe can be
!      used in a simulation for a region 30N to 50N and 130W to 70W -- the 
!      code will extract the appropriate data. The two grids do not have to
!      be the same resolution. Area-averaging will work for full => partial
!      grid but obviously will not work for partial => full grid.
!
! o 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 as
!
!   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)
!
! o Error checks:
!
!   Overlap weights of input cells sum to 1 for each output cell.
!   Global sum of dummy field is conserved for input => output area-average.
! 
! 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 longitude points
    integer , intent(in)    :: nlat_i                 !input  grid: number of latitude  points
    integer , intent(in)    :: numlon_i(nlat_i)       !input  grid: number lon points at each lat
    real(r8), intent(in)    :: lon_i(nlon_i+1,nlat_i) !input grid: longitude, west edge (degrees)
    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: latitude, south edge (degrees)
    real(r8), intent(in)    :: area_i(nlon_i,nlat_i)  !input grid: cell area
    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 lon points at each lat
    real(r8), intent(in)    :: lon_o(nlon_o+1,nlat_o) !output grid: longitude, west edge  (degrees)
    real(r8), intent(in)    :: lat_o(nlat_o+1)        !output grid: latitude, south edge (degrees)
    real(r8), intent(in)    :: area_o                 !output grid: cell area
    real(r8), intent(in)    :: fland_o                !output grid: fraction that is land
    integer , intent(out)   :: novr_i2o               !number of overlapping input cells
    integer , intent(out)   :: iovr_i2o(maxovr)        !lon index of overlap input cell
    integer , intent(out)   :: jovr_i2o(maxovr)        !lat index of overlap input cell
    real(r8), intent(out)   :: wovr_i2o(maxovr)        !weight    of overlap input cell
! --------------------------------------------------------------------

! ------------------------ local variables ---------------------------
    real(r8) :: relerr = 0.00001       !relative error for error checks
    integer  :: ii                     !input  grid longitude loop index
    integer  :: ji                     !input  grid latitude  loop index
    integer  :: n                      !overlap index
    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
! --------------------------------------------------------------------

! -----------------------------------------------------------------